Line data Source code
1 : !!****m* ABINIT/m_dmft
2 : !! NAME
3 : !! m_dmft
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 : !! INPUTS
14 : !!
15 : !! OUTPUT
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : ! nvtx related macro definition
26 : #include "nvtx_macros.h"
27 :
28 : MODULE m_dmft
29 :
30 : use defs_abitypes
31 : use defs_basis
32 : !use netcdf
33 : use m_xmpi
34 : use m_abicore
35 : use m_data4entropyDMFT
36 : use m_errors
37 :
38 : use m_crystal, only : crystal_t
39 : use m_datafordmft, only : chipsi_print,chipsi_renormalization,compute_wannier,print_wannier
40 : use m_dftu_self, only : dftu_self
41 : use m_energy, only : compute_dftu_energy,compute_energy,compute_free_energy,&
42 : & destroy_energy,energy_type,init_energy
43 : use m_forctqmc, only : ctqmc_calltriqs_c,qmc_prep_ctqmc
44 : use m_green, only : check_fourier_green,compute_green,copy_green,destroy_green,destroy_green_tau, &
45 : & fermi_green,fourier_green,green_type,icip_green,init_green,init_green_tau,integrate_green, &
46 : & local_ks_green,print_green,printocc_green
47 : use m_hu, only : destroy_hu,hu_type,init_hu
48 : use m_hubbard_one, only : hubbard_one
49 : use m_matlu, only : add_matlu,copy_matlu,destroy_matlu,diff_matlu,identity_matlu,init_matlu,inverse_matlu, &
50 : & matlu_type,print_matlu,sym_matlu,xmpi_matlu
51 : use m_oper, only : destroy_oper,diff_oper,downfold_oper,gather_oper,init_oper,inverse_oper,oper_type,trace_oper
52 : use m_paw_dmft, only : paw_dmft_type
53 : use m_pawang, only : pawang_type
54 : use m_pawtab, only : pawtab_type
55 : use m_self, only : dc_self,destroy_self,initialize_self,new_self,print_self,rw_self,self_type
56 : use m_time, only : timab
57 :
58 : #ifdef HAVE_GPU_MARKERS
59 : use m_nvtx_data
60 : #endif
61 :
62 : implicit none
63 :
64 : private
65 :
66 : public :: dmft_solve
67 : public :: impurity_solve
68 : public :: dyson
69 : public :: spectral_function
70 : !!***
71 :
72 : contains
73 :
74 : !!****f* ABINIT/dmft_solve
75 : !! NAME
76 : !! dmft_solve
77 : !!
78 : !! FUNCTION
79 : !! Solve the DMFT loop from PAW data.
80 : !!
81 : !! INPUTS
82 : !! cryst_struc <type(crystal_t)>=crystal structure data
83 : !! istep = iteration step of the DFT+DMFT self-consistent cycle.
84 : !! dft_occup <type(oper_type)> = DFT occupations numbers of the correlated orbitals
85 : !! mpi_enreg=information about MPI parallelization
86 : !! paw_dmft <type(paw_dmft_type)> = data for self-consistent DFT+DMFT calculations.
87 : !! pawang <type(pawang)>=paw angular mesh and related data
88 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
89 : !! pawprtvol = option for printing
90 : !!
91 : !! OUTPUT
92 : !! paw_dmft <type(paw_dmft_type)> = data for self-consistent DFT+DMFT calculations.
93 : !!
94 : !! NOTES
95 : !!
96 : !! SOURCE
97 :
98 93 : subroutine dmft_solve(cryst_struc,istep,dft_occup,mpi_enreg,paw_dmft,pawang,pawtab,pawprtvol)
99 :
100 : !Arguments ------------------------------------
101 : integer, intent(in) :: istep,pawprtvol
102 : type(MPI_type), intent(in) :: mpi_enreg
103 : type(pawang_type), intent(in) :: pawang
104 : type(crystal_t), intent(in) :: cryst_struc
105 : type(paw_dmft_type), intent(inout) :: paw_dmft
106 : type(pawtab_type), intent(inout) :: pawtab(paw_dmft%ntypat)
107 : type(oper_type), intent(in) :: dft_occup
108 : !Local variables ------------------------------
109 : integer :: check,dmft_iter,idmftloop,istep_iter,itypat,myproc,natom
110 : integer :: ntypat,opt_diff,opt_maxent,opt_moments,opt_renorm,prtopt
111 : !logical :: etot_var
112 : logical :: dmft_optim,t2g,x2my2d
113 : real(dp) :: tsec(2)
114 : character(len=200) :: char_enddmft
115 93 : type(green_type) :: green,greendft,weiss
116 93 : type(self_type) :: self,self_new
117 93 : type(energy_type) :: energies_dmft,energies_tmp
118 93 : type(oper_type) :: identity_oper,oper_tmp
119 93 : type(hu_type), allocatable :: hu(:)
120 : character(len=4) :: part2,part3
121 : !character(len=5) :: thdyn
122 : character(len=500) :: message
123 : !************************************************************************
124 :
125 : DBG_ENTER('COLL')
126 : ABI_NVTX_START_RANGE(NVTX_DMFT_SOLVE)
127 :
128 93 : myproc = paw_dmft%myproc
129 93 : check = paw_dmft%dmftcheck ! checks enabled
130 93 : t2g = (paw_dmft%dmft_t2g == 1)
131 93 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
132 93 : natom = paw_dmft%natom
133 93 : ntypat = paw_dmft%ntypat
134 93 : dmft_iter = paw_dmft%dmft_iter
135 93 : opt_maxent = paw_dmft%dmft_prt_maxent
136 93 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) ! skip some unneeded calls to compute_green with TRIQS
137 : !paw_dmft%dmft_fermi_prec=tol5
138 : !paw_dmft%dmft_fermi_prec = paw_dmft%dmft_charge_prec * ten
139 : !paw_dmft%dmft_charge_prec=20_dp ! total number of electron.
140 : !paw_dmft%dmft_prgn=1
141 93 : paw_dmft%dmft_prgn = 0
142 : !etot_var = .true.
143 : !thdyn="fcalc"
144 : !thdyn = "ecalc"
145 : !if (thdyn == "ecalc") then ! valid
146 93 : part2 = "both"
147 93 : part3 = "none"
148 : !else if (thdyn == "fcalc") then ! not tested
149 : ! part2 = "corr"
150 : ! part3 = "band"
151 : !end if
152 :
153 93 : opt_moments = merge(1,0,paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7)
154 93 : prtopt = merge(2,0,dmft_optim)
155 93 : opt_diff = merge(1,0,dmft_optim)
156 :
157 93 : if (check == 1) then
158 1 : write(message,'(2a)') ch10,' DMFT Checks are enabled '
159 : else
160 92 : write(message,'(2a)') ch10,' DMFT Checks will not be performed'
161 : end if ! check
162 93 : call wrtout(std_out,message,'COLL')
163 :
164 93 : if (istep == 0) then
165 0 : message = ' istep should not be equal to zero'
166 0 : ABI_BUG(message)
167 : end if
168 :
169 : !spaceComm=paw_dmft%spacecomm
170 : !if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
171 : !call xmpi_barrier(spaceComm)
172 :
173 93 : call initialize_self(self,paw_dmft,opt_moments=opt_moments)
174 93 : call init_energy(energies_dmft,natom)
175 :
176 : !===========================================================================
177 : !== First construct DFT green function (Init, Compute, Integrate, Print)
178 : !===========================================================================
179 93 : write(message,'(6a)') ch10," ==========================================================================", &
180 93 : & ch10," ===== Check: DFT Green's Function Calculation with unnormalized orbitals",&
181 186 : & ch10," =========================================================================="
182 93 : call wrtout(std_out,message,'COLL')
183 93 : call icip_green("DFT",greendft,paw_dmft,3,self,opt_moments=opt_moments)
184 : !call print_green('DFT_NOT_renormalized',greendft,1,paw_dmft,pawprtvol=1,opt_wt=1)
185 :
186 : !== Compare greendft%occup and dft_occup: check that DFT green function is fine
187 : !----------------------------------------------------------------------
188 93 : write(message,'(2a)') ch10," == Compare local occupations from DFT Green's function &
189 186 : &with the downfold of the Fermi-Dirac occupations =="
190 93 : call wrtout(std_out,message,'COLL')
191 :
192 93 : if(paw_dmft%dmft_magnfield .gt. 0) then
193 13 : write(message, '(2a,a)') ch10, 'Warning: Check in local occupation is removed due to applied magnetic field'
194 13 : call wrtout(std_out,message,'COLL')
195 : else
196 : call diff_oper("occupations from DFT Green's function","Fermi-Dirac occupations", &
197 80 : & greendft%occup,dft_occup,1,paw_dmft%dmft_tolfreq)
198 : endif
199 : ! write(message,'(2a)') ch10,&
200 : !& ' ***** => Warning : diff_oper is suppressed for test'
201 : ! call wrtout(std_out,message,'COLL')
202 93 : write(message,'(2a)') ch10," ***** => Calculation of DFT Green's function is thus correct ****"
203 93 : call wrtout(std_out,message,'COLL')
204 93 : call destroy_green(greendft)
205 :
206 : !== Orthonormalize chipsi
207 : !----------------------------------------------------------------------
208 93 : call timab(621,1,tsec(:))
209 : !natomcor=0
210 : !do iatom=1,paw_dmft%natom
211 : ! if(paw_dmft%lpawu(iatom).ne.-1) then
212 : ! natomcor=natomcor+1
213 : ! end if
214 : !end do
215 : opt_renorm = merge(2,paw_dmft%dmft_wanorthnorm,paw_dmft%nspinor==2.and. &
216 93 : &(paw_dmft%dmft_solv == 9 ))
217 : ! &(paw_dmft%dmft_solv == 8 .or. paw_dmft%dmft_solv == 9 )) !.or. paw_dmft%dmft_solv == 10))
218 :
219 93 : if (paw_dmft%dmft_solv /= -1) then
220 92 : call chipsi_renormalization(paw_dmft,opt=opt_renorm)
221 92 : if (paw_dmft%dmft_prtwan == 1) then
222 0 : call compute_wannier(paw_dmft,mpi_enreg)
223 0 : if (myproc == 0) then
224 0 : call print_wannier(paw_dmft,istep)
225 : end if
226 0 : ABI_FREE(paw_dmft%wannier)
227 : end if ! dmft_prtwan
228 :
229 92 : write(message,'(2a)') ch10,' == Check downfold(upfold)=identity =='
230 92 : call wrtout(std_out,message,'COLL')
231 :
232 : ! Check that downfold_oper(upfold_oper)=I
233 92 : call init_oper(paw_dmft,identity_oper,opt_ksloc=2)
234 92 : call init_oper(paw_dmft,oper_tmp,opt_ksloc=2)
235 92 : call identity_matlu(identity_oper%matlu(:),natom)
236 92 : call downfold_oper(oper_tmp,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=4)
237 92 : call xmpi_matlu(oper_tmp%matlu(:),natom,paw_dmft%distrib%comm_kpt)
238 92 : call sym_matlu(oper_tmp%matlu(:),paw_dmft)
239 92 : call diff_matlu("Downfold(Upfold)","Identity",oper_tmp%matlu(:),identity_oper%matlu(:),natom,0,tol4)
240 92 : call destroy_oper(oper_tmp)
241 92 : call destroy_oper(identity_oper)
242 :
243 : ! ===========================================================================
244 : ! == re-construct DFT green function with new chipsis
245 : ! ===========================================================================
246 : write(message,'(6a)') &
247 92 : & ch10," ========================================================================", &
248 92 : & ch10," ===== Check: DFT Green's Function Calculation with normalized orbitals", &
249 184 : & ch10," ========================================================================"
250 : end if ! dmft_solv/=1
251 93 : call timab(621,2,tsec(:))
252 93 : call wrtout(std_out,message,'COLL')
253 :
254 93 : call icip_green("DFT renormalized",greendft,paw_dmft,pawprtvol,self,opt_moments=opt_moments,opt_log=paw_dmft%dmft_triqs_entropy)
255 : !call print_green('DFT_renormalized',greendft,1,paw_dmft,pawprtvol=1,opt_wt=1)
256 :
257 : !== Define Interaction from input upawu and jpawu
258 : !----------------------------------------------------------------------
259 482 : ABI_MALLOC(hu,(ntypat))
260 93 : call init_hu(hu(:),paw_dmft,pawtab(:))
261 :
262 93 : call dc_self(greendft%charge_matlu(:,:),self%hdc%matlu(:),hu(:),paw_dmft,pawtab(:),greendft%occup%matlu(:))
263 :
264 : ! Need to store idmftloop and set it to zero to avoid useless print_energy in ab_out
265 93 : idmftloop = paw_dmft%idmftloop
266 93 : paw_dmft%idmftloop = 0
267 93 : call compute_energy(energies_dmft,greendft,paw_dmft,pawprtvol,pawtab(:),self,occ_type=" lda",part='both')
268 93 : if (paw_dmft%dmft_triqs_entropy == 1) then
269 0 : call compute_free_energy(energies_dmft,paw_dmft,greendft,"band")
270 : end if
271 93 : paw_dmft%idmftloop = idmftloop
272 :
273 93 : if ((paw_dmft%dmft_prgn == 1) .and. (paw_dmft%lchipsiortho == 1)) then
274 0 : call local_ks_green(greendft,paw_dmft,prtopt=1)
275 : end if
276 : !call printocc_green(greendft,9,paw_dmft,3,chtype="DFT GREEN PSICHI")
277 :
278 : write(message,'(7a)') &
279 93 : & ch10,' =============================', &
280 93 : & ch10,' ===== Define self-energy', &
281 186 : & ch10,' =============================',ch10
282 93 : call wrtout(std_out,message,'COLL')
283 :
284 : ! Set Hu in density representation for calculation of entropy if needed...
285 93 : if (paw_dmft%dmft_entropy > 0) then
286 30 : do itypat=1,ntypat
287 21 : if (hu(itypat)%lpawu == -1) cycle
288 912 : call data4entropyDMFT_setHu(paw_dmft%forentropyDMFT,itypat,dble(hu(itypat)%udens(:,:)))
289 : end do ! itypat
290 : end if ! dmft_entropy=1
291 :
292 : !== define self from scratch or file and double counting
293 : !----------------------------------------------------------------------
294 : !- Self allocated
295 :
296 : !- Read self or do self=hdc
297 : !if(paw_dmft%dmft_solv==4) then
298 : ! write(std_out,*) "shift before rw_self",self%qmc_shift(1)
299 : ! call make_qmcshift_self(cryst_struc,hu,self)
300 : !end if
301 93 : call timab(627,1,tsec(:))
302 93 : call rw_self(self,paw_dmft,2,opt_rw=1,istep_iter=1000*istep)
303 93 : call timab(627,2,tsec(:))
304 :
305 : !== If QMC is used, and self energy is read for file, then
306 : !== one does NOT shifts the self-energy because it was already shifted when writed,
307 : !== and thus then weiss will be shifted
308 : !----------------------------------------------------------------------
309 : !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf==1) &
310 : !& call make_qmcshift_self(cryst_struc,hu,self)
311 : !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf/=1) &
312 : !& call make_qmcshift_self(cryst_struc,hu,self,apply=.true.)
313 :
314 93 : call destroy_green(greendft) ! destroy DFT green function
315 93 : call print_self(self,"print_dc",paw_dmft,2)
316 :
317 : !===========================================================================
318 : !== Construct green function with the self-energy.
319 : !===========================================================================
320 : write(message,'(6a)') &
321 93 : & ch10," ===================================================================", &
322 93 : & ch10," ===== Green's Function Calculation with input self-energy ========", &
323 186 : & ch10," ==================================================================="
324 93 : call wrtout(std_out,message,'COLL')
325 93 : if (dmft_optim) then
326 0 : call init_green(green,paw_dmft,opt_moments=opt_moments)
327 : else
328 93 : call icip_green("DFT+DMFT",green,paw_dmft,pawprtvol,self,opt_self=1,opt_moments=opt_moments)
329 : !call print_green('beforefermi_green',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
330 : ! call abi_abort('COLL')
331 : end if
332 :
333 : !== Find fermi level
334 : !---------------------------------------------------------------------
335 : !write(message,'(2a,i3,13x,a)') ch10,' === Compute green function from self-energy'
336 :
337 93 : call fermi_green(green,paw_dmft,self)
338 93 : call compute_green(green,paw_dmft,0,self,opt_self=1,opt_nonxsum=1,opt_restart_moments=1)
339 93 : call integrate_green(green,paw_dmft,prtopt)
340 :
341 93 : if (dmft_optim) then
342 0 : call printocc_green(green,5,paw_dmft,3,chtype="DFT+DMFT")
343 : end if
344 :
345 : !== define weiss field only for the local quantities (opt_oper=2)
346 : !----------------------------------------------------------------------
347 : ! write(std_out,*) "nkpt befreo init_greenweiss",ifreq,paw_dmft%nkpt
348 93 : call init_green(weiss,paw_dmft,opt_oper_ksloc=2,opt_moments=opt_moments,opt_moments_ksloc=2,opt_occup_ksloc=2)
349 : ! do ifreq=1,weiss%nw
350 : ! write(std_out,*) "nkpt from weiss1",ifreq,weiss%oper(ifreq)%nkpt
351 : ! enddo
352 :
353 : !== Check fourier transforms
354 : !----------------------------------------------------------------------
355 93 : if (check == 1) then
356 1 : call check_fourier_green(cryst_struc,green,paw_dmft)
357 : end if
358 :
359 : !== If QMC is used, and self energy is not read for file, then
360 : !== one shifts the self-energy, and thus then weiss will be shifted
361 : !== after dyson, in a coherent way regarding qmc_shift and qmc_xmu.
362 : !----------------------------------------------------------------------
363 : !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf/=1) &
364 : !& call make_qmcshift_self(cryst_struc,hu,self,apply=.true.)
365 : !if(paw_dmft%dmft_solv==4) write(std_out,*) "shift after make_qmcshift_self",self%qmc_shift(1)
366 :
367 : write(message,'(6a)') &
368 93 : & ch10,' ======================================================', &
369 93 : & ch10,' ===== DMFT Loop starts here ========', &
370 186 : & ch10,' ======================================================'
371 93 : call wrtout(std_out,message,'COLL')
372 :
373 : ABI_NVTX_START_RANGE(NVTX_DMFT_SOLVE_LOOP)
374 : !=======================================================================
375 : !=== dmft loop =======================================================
376 199 : do idmftloop=1,dmft_iter
377 : !paw_dmft%idmftloop=idmftloop
378 106 : paw_dmft%idmftloop = paw_dmft%idmftloop + 1
379 : ! =======================================================================
380 106 : istep_iter = 1000*istep + idmftloop
381 :
382 106 : write(message,'(2a,i3,13x,a)') ch10,&
383 212 : & ' ===== DMFT Loop : ITER number',paw_dmft%idmftloop,'========'
384 106 : call wrtout(std_out,message,'COLL')
385 :
386 : ! == Dyson Equation G,self -> weiss(w)
387 : ! ---------------------------------------------------------------------
388 106 : call dyson(green,paw_dmft,self,weiss,opt_weissself=1)
389 : ! call print_green('afterDyson',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
390 : ! call abi_abort('COLL')
391 :
392 : ! == Printout local "occupations" from weiss field (useless)
393 106 : if (abs(pawprtvol) > 3 .and. opt_moments == 0) then
394 0 : call integrate_green(weiss,paw_dmft,2,opt_ksloc=2)
395 0 : call printocc_green(weiss,5,paw_dmft,3,opt_weissgreen=1)
396 : end if
397 :
398 : ! === Prepare data, solve Impurity problem: weiss(w) -> G(w)
399 : ! ---------------------------------------------------------------------
400 106 : call initialize_self(self_new,paw_dmft,opt_moments=opt_moments)
401 :
402 106 : call impurity_solve(cryst_struc,green,hu(:),paw_dmft,pawang,pawtab(:),self,self_new,weiss,pawprtvol) ! weiss-> green, or self if dmft_solv=1
403 : ! if(paw_dmft%dmft_solv==4) write(std_out,*) "shift after impurity",self%qmc_shift(1)
404 :
405 : ! == Compute double counting from charge from green_solver
406 : ! ---------------------------------------------------------------------
407 130 : if (green%has_charge_matlu_solver /= 2) green%charge_matlu_solver(:,:) = green%charge_matlu(:,:)
408 :
409 106 : if (paw_dmft%dmft_solv >= 5) then
410 34 : call dc_self(green%charge_matlu_solver(:,:),self_new%hdc%matlu(:),hu(:),paw_dmft,pawtab(:),green%occup_tau%matlu(:))
411 : else
412 72 : call dc_self(green%charge_matlu_solver(:,:),self_new%hdc%matlu(:),hu(:),paw_dmft,pawtab(:),green%occup%matlu(:))
413 : end if
414 :
415 106 : if (abs(paw_dmft%dmft_solv) >= 5) then
416 34 : call destroy_green_tau(green)
417 : end if
418 :
419 : ! == Solve dyson equation. G_imp(w), weiss_imp(w) -> Self_imp(w)
420 : ! ---------------------------------------------------------------------
421 : ! if dmft_solv==1, self is computed previously
422 106 : if (abs(paw_dmft%dmft_solv) /= 1) then
423 104 : call dyson(green,paw_dmft,self_new,weiss,opt_weissself=2)
424 : end if
425 : ! do ifreq=1,green%nw
426 : ! call sym_matlu(cryst_struc,self%oper(ifreq)%matlu,pawang)
427 : ! enddo
428 :
429 : ! == Possibility if imposing self (opt_rw==3)
430 : ! ---------------------------------------------------------------------
431 106 : call timab(627,1,tsec(:))
432 106 : call rw_self(self_new,paw_dmft,prtopt=2,opt_rw=3,istep_iter=istep_iter)
433 106 : call timab(627,2,tsec(:))
434 :
435 : ! Print dc computed just before and self computed before in dyson or
436 : ! impurity_solve
437 106 : if (abs(pawprtvol) >= 3) then
438 70 : write(message,'(2a)') ch10," == Old self (before impurity solver)"
439 70 : call wrtout(std_out,message,'COLL')
440 70 : call print_self(self,"print_dc",paw_dmft,2)
441 70 : write(message,'(2a)') ch10," == New self (from impurity solver)"
442 70 : call wrtout(std_out,message,'COLL')
443 70 : call print_self(self_new,"print_dc",paw_dmft,2)
444 : end if ! abs(pawprtvol)>=3
445 :
446 : ! if(paw_dmft%dmft_solv==4) write(std_out,*) "shift before computeenergy ",self%qmc_shift(1)
447 : ! == Compute Energy with NEW self-energy and edc from green_solver,
448 : ! new local green function and old occupations for eband
449 : ! fermi level not optimized for this self_energy.
450 : ! ---------------------------------------------------------------------
451 : ! green= local green function and local charge comes directly from solver
452 : ! green= ks green function and occupations comes from old_self
453 106 : call compute_energy(energies_dmft,green,paw_dmft,pawprtvol,pawtab(:),self_new,occ_type="nlda",part=part2)
454 106 : if (paw_dmft%dmft_triqs_entropy == 1) then
455 0 : call compute_free_energy(energies_dmft,paw_dmft,green,"impu",self_new)
456 : end if
457 :
458 : ! == Mix new and old self_energies and double countings
459 : ! ---------------------------------------------------------------------
460 106 : write(message,'(3a)') ch10," == Linear mixing of old and new self-energy and double counting",ch10
461 106 : call wrtout(std_out,message,'COLL')
462 106 : call new_self(self,self_new,paw_dmft) ! self,self_new => self
463 106 : write(message,'(2a)') ch10," == After mixing,"
464 : !print *, " my_rank newself", my_rank,self%oper(1)%matlu(1)%mat(1,1,1,1,1)
465 106 : call wrtout(std_out,message,'COLL')
466 106 : call print_self(self,"print_dc",paw_dmft,2) ! print self and DC
467 106 : call destroy_self(self_new)
468 :
469 : ! == Compute green function self -> G(k)
470 : ! ---------------------------------------------------------------------
471 106 : if (.not. dmft_optim) then
472 106 : call compute_green(green,paw_dmft,1,self,opt_self=1,opt_nonxsum=1)
473 106 : call integrate_green(green,paw_dmft,3,opt_diff=1) !,opt_nonxsum=1)
474 :
475 106 : call printocc_green(green,5,paw_dmft,3,chtype="DFT+DMFT")
476 : ! call printocc_green(green,9,paw_dmft,3,chtype="DMFT FULL")
477 106 : if(paw_dmft%lchipsiortho == 1 .and. paw_dmft%dmft_prgn == 1) then
478 0 : call local_ks_green(green,paw_dmft,prtopt=1)
479 : end if
480 : end if ! dmft_optim=0
481 :
482 : ! == Find fermi level
483 : ! ---------------------------------------------------------------------
484 106 : call fermi_green(green,paw_dmft,self)
485 106 : call compute_green(green,paw_dmft,0,self,opt_self=1,opt_nonxsum=1,opt_log=paw_dmft%dmft_triqs_entropy,opt_restart_moments=1)
486 106 : call integrate_green(green,paw_dmft,prtopt,opt_diff=opt_diff,opt_ksloc=3,opt_fill_occnd=1)
487 :
488 106 : if (dmft_optim) then
489 0 : call printocc_green(green,5,paw_dmft,3,chtype="DFT+DMFT")
490 : end if
491 :
492 : ! call abi_abort('COLL')
493 :
494 : ! == Compute Energy with Mixed self-energy and green function recomputed with new self
495 : ! ---------------------------------------------------------------------
496 : ! green= lattice green function computed from self for a given chemical potential mu (self_mixed,mu)
497 : ! green= local green function is computed from lattice green function(self_mixed,mu)
498 : ! green= occupations are computed with lattice green function(self_mixed,mu)
499 106 : call compute_energy(energies_dmft,green,paw_dmft,pawprtvol,pawtab(:),self,occ_type="nlda",part=part3)
500 :
501 : ! == Save self on disk
502 : ! ---------------------------------------------------------------------
503 106 : call timab(627,1,tsec(:))
504 106 : call rw_self(self,paw_dmft,prtopt=2,opt_rw=2,opt_maxent=opt_maxent)
505 106 : call timab(627,2,tsec(:))
506 :
507 : ! == Test convergency
508 : ! ---------------------------------------------------------------------
509 106 : char_enddmft = "DFT+DMFT (end of DMFT loop)"
510 199 : if (green%ifermie_cv == 1 .and. self%iself_cv == 1 .and. green%ichargeloc_cv == 1 .and. paw_dmft%idmftloop > 1) then
511 0 : write(message,'(a,8x,a)') ch10,"DMFT Loop is converged !"
512 0 : call wrtout(std_out,message,'COLL')
513 0 : char_enddmft = "converged DMFT"
514 0 : exit
515 : end if
516 : ! =======================================================================
517 : ! === end dmft loop ====================================================
518 : end do ! idmftloop
519 : ABI_NVTX_END_RANGE()
520 : !=========================================================================
521 :
522 : !== Save self on disk
523 : !-------------------------------------------------------------------------
524 93 : if (.not. dmft_optim) then
525 93 : call timab(627,1,tsec(:))
526 93 : call rw_self(self,paw_dmft,prtopt=2,opt_rw=2)
527 93 : call timab(627,2,tsec(:))
528 : end if
529 :
530 : !paw_dmft%idmftloop=0
531 :
532 93 : write(message,'(2a,13x,a)') ch10,' ===== DMFT Loop : END ','========'
533 93 : call wrtout(std_out,message,'COLL')
534 :
535 93 : if (paw_dmft%dmft_entropy >= 1) then
536 : ! compute Edc for U=1 and J=U/J
537 9 : call init_energy(energies_tmp,natom)
538 : !call compute_dftu_energy(cryst_struc,energies_tmp,green,paw_dmft,pawtab)
539 9 : call compute_dftu_energy(energies_tmp,green,paw_dmft,pawtab(:),paw_dmft%forentropyDMFT%J_over_U)
540 9 : call data4entropyDMFT_setDc(paw_dmft%forentropyDMFT,energies_tmp%e_dc(:))
541 9 : call destroy_energy(energies_tmp,paw_dmft)
542 : end if ! dmft_entropy=1
543 :
544 : !== Compute final values for green functions, occupations, and spectral function
545 : !--------------------------------------------------------------------------------
546 : !Do not compute here, because, one want a energy computed after the
547 : !solver (for Hubbard I and DFT+U).
548 93 : if (.not. dmft_optim) then
549 93 : call compute_green(green,paw_dmft,1,self,opt_self=1,opt_nonxsum=1)
550 93 : call integrate_green(green,paw_dmft,2,opt_fill_occnd=1) !,opt_nonxsum=1)
551 : end if
552 : !call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,opt=0)
553 93 : idmftloop = paw_dmft%idmftloop
554 93 : paw_dmft%idmftloop = 0
555 93 : call compute_energy(energies_dmft,green,paw_dmft,pawprtvol,pawtab(:),self,occ_type="nlda",part="band")
556 93 : paw_dmft%idmftloop = idmftloop
557 93 : if (paw_dmft%dmft_triqs_entropy == 1) then
558 0 : call compute_free_energy(energies_dmft,paw_dmft,green,"main",self)
559 : end if
560 :
561 : !write(message,'(2a,13x,a)') ch10,' ===== DMFT Loop is finished'
562 : !call wrtout(ab_out,message,'COLL')
563 : !write(std_out,*) "PRINTOCC INITIAL"
564 93 : call printocc_green(green,9,paw_dmft,3,chtype=char_enddmft)
565 : !write(std_out,*) "KS=czero"
566 : !green%occup%ks=czero
567 : !write(std_out,*) "PRINTOCC AFTER KS=0"
568 : !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT")
569 : !write(std_out,*) "UPFOLD_OPER"
570 : !call upfold_oper(green%occup,paw_dmft,1)
571 : !write(std_out,*) "PRINTOCC AFTER UPFOLD_OPER"
572 : !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT")
573 : !write(std_out,*) "MATLU=czero"
574 : !green%occup%matlu(1)%mat=czero
575 : !green%occup%ks(:,:,:,:)=cmplx(real(green%occup%ks(:,:,:,:)))
576 : !write(std_out,*) "PRINTOCC AFTER MATLU=0 AND IMAG KS=0"
577 : !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT")
578 : !write(std_out,*) "LOC_OPER"
579 : !call loc_oper(green%occup,paw_dmft,1)
580 : !write(std_out,*) "PRINTOCC AFTER LOC_OPER"
581 : !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT")
582 : !call flush_unit(std_out)
583 : !call abi_abort('COLL')
584 93 : if (paw_dmft%dmft_solv <= 2 .and. paw_dmft%prtdos >= 1) then
585 0 : call spectral_function(cryst_struc,green,hu(:),paw_dmft,pawtab(:),self,pawprtvol)
586 : end if
587 93 : call destroy_green(weiss)
588 93 : call destroy_green(green)
589 : !todo_ab rotate back density matrix into unnormalized basis just for
590 : !printout
591 93 : call destroy_hu(hu(:),ntypat)
592 93 : call destroy_self(self)
593 93 : call destroy_energy(energies_dmft,paw_dmft)
594 :
595 93 : write(message,'(2a,13x,a)') ch10,' ===== DMFT : END ','========'
596 93 : call wrtout(std_out,message,'COLL')
597 :
598 296 : ABI_FREE(hu)
599 :
600 : ABI_NVTX_END_RANGE()
601 : DBG_EXIT("COLL")
602 :
603 93 : end subroutine dmft_solve
604 : !!***
605 :
606 : !!****f* ABINIT/impurity_solve
607 : !! NAME
608 : !! impurity_solve
609 : !!
610 : !! FUNCTION
611 : !! Solve the Impurity problem
612 : !!
613 : !! INPUTS
614 : !! cryst_struc <type(crystal_t)>=crystal structure data
615 : !! green <type(green_type)>= green function data
616 : !! hu <type(hu_type)>= U interaction
617 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
618 : !! pawang <type(pawang)>=paw angular mesh and related data
619 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
620 : !! self_old,self_new <type(self_type)>= variables related to self-energy
621 : !! weiss <type(green_type)>= weiss function data
622 : !! pawprtvol = option for printing
623 : !!
624 : !! OUTPUT
625 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
626 : !!
627 : !! NOTES
628 : !!
629 : !! SOURCE
630 :
631 106 : subroutine impurity_solve(cryst_struc,green,hu,paw_dmft,pawang,pawtab,&
632 : & self_old,self_new,weiss,pawprtvol)
633 :
634 : !Arguments ------------------------------------
635 : !scalars
636 : ! type(pawang_type), intent(in) :: pawang
637 : type(crystal_t), intent(in) :: cryst_struc
638 : type(green_type), intent(inout) :: green,weiss
639 : type(paw_dmft_type), intent(inout) :: paw_dmft
640 : type(hu_type), intent(inout) :: hu(paw_dmft%ntypat)
641 : !type(MPI_type), intent(in) :: mpi_enreg
642 : type(pawang_type), intent(in) :: pawang
643 : type(pawtab_type), intent(in) :: pawtab(paw_dmft%ntypat)
644 : type(self_type), intent(inout) :: self_new,self_old
645 : integer, intent(in) :: pawprtvol
646 : !Local variables ------------------------------
647 : real(dp) :: tsec(2)
648 : character(len=500) :: message
649 : ! integer iatom,il,i_nd,isppol,lpawu,im,Nd,nrat,nsweeptot
650 : ! real(dp) :: acc,kx
651 : ! real(dp), allocatable :: correl(:,:),g0(:,:),gtmp(:,:)
652 : !scalars
653 : !************************************************************************
654 : !character(len=500) :: message
655 :
656 106 : call timab(622,1,tsec(:))
657 : ABI_NVTX_START_RANGE(NVTX_DMFT_IMPURITY_SOLVE)
658 : !=======================================================================
659 : !== Prepare data for Hirsch Fye QMC
660 : !== NB: for CTQMC, Fourier Transformation are done inside the CTQMC code
661 : !=======================================================================
662 : !if(abs(paw_dmft%dmft_solv)==4) then
663 : ! == Initialize weiss and green functions for fourier transformation
664 : ! -------------------------------------------------------------------
665 : ! write(message,'(2a,i3,13x,a)') ch10,' === Initialize Weiss field G_0(tau)'
666 : ! call wrtout(std_out,message,'COLL')
667 : ! call init_green_tau(weiss,paw_dmft)
668 : ! call init_green_tau(green,paw_dmft)
669 : ! in init_solver
670 :
671 : ! == Print weiss function G_0(tau=0-) before computation (really useless check)
672 : ! ------------------------------------------------------------------------------
673 : ! if(abs(pawprtvol)>3) then
674 : ! write(message,'(2a,i3,13x,a)') ch10,' === Check G_0(tau=0-) first'
675 : ! call wrtout(std_out,message,'COLL')
676 : ! call printocc_green(weiss,6,paw_dmft,3)
677 : ! end if
678 :
679 : ! == Fourier transform of weiss Field
680 : ! ------------------------------------
681 : ! for fourier of KS green functions
682 : ! call fourier_green(cryst_struc,weiss,mpi_enreg,paw_dmft,pawang,pawtab,1)
683 : ! write(message,'(2a,i3,13x,a)') ch10,' === Inverse Fourier Transform w->t of Weiss Field'
684 : ! call wrtout(std_out,message,'COLL')
685 : ! call fourier_green(cryst_struc,weiss,paw_dmft,pawang,opt_ksloc=2,opt_tw=-1)
686 :
687 : ! == Print weiss function G2_0(tau=0-)
688 : ! --------------------------------------
689 : ! call printocc_green(weiss,6,paw_dmft,3,opt_weissgreen=1)
690 :
691 : ! for fourier of KS green functions
692 : ! call fourier_green(cryst_struc,weiss,mpi_enreg,paw_dmft,pawang,pawtab,1)
693 : ! == Print G_0(tau) in files
694 : ! ---------------------------
695 : ! if(paw_dmft%dmft_prgn==1) then
696 : ! call print_green('weiss',weiss,1,paw_dmft,pawprtvol=1,opt_wt=2)
697 : ! end if
698 :
699 106 : if (abs(paw_dmft%dmft_solv) >= 5) then
700 : ! == Initialize green functions for imaginary times
701 : ! -------------------------------------------------------------------
702 34 : write(message,'(2a)') ch10," === Initialize Green's function G(tau)"
703 34 : call wrtout(std_out,message,'COLL')
704 34 : call init_green_tau(green,paw_dmft)
705 :
706 : end if
707 : !=======================================================================
708 : !== End preparation of QMC
709 : !=======================================================================
710 :
711 : !=======================================================================
712 : !== Solve impurity model =============================================
713 : !=======================================================================
714 106 : write(message,'(2a)') ch10,' === Solving impurity model'
715 106 : call wrtout(std_out,message,'COLL')
716 106 : if (abs(paw_dmft%dmft_solv) == 1) then
717 :
718 : ! == DFT+U for test -> self
719 : ! -------------------
720 2 : call dftu_self(cryst_struc,green,paw_dmft,pawtab(:),self_new,opt_dftu=1,prtopt=pawprtvol)
721 :
722 104 : else if (abs(paw_dmft%dmft_solv) == 2) then
723 :
724 : ! == Hubbard One -> green
725 : ! -------------------
726 69 : call hubbard_one(cryst_struc,green,hu(:),paw_dmft,pawprtvol,self_old%hdc,weiss)
727 :
728 : !else if(abs(paw_dmft%dmft_solv)==4) then
729 :
730 : ! == QMC
731 : ! -------------------
732 : ! call copy_green(weiss,green,opt_tw=1)
733 : ! call qmc_prep
734 : ! message = ' === QMC not yet distributed '
735 : ! ABI_ERROR(message)
736 : ! call qmc_prep(cryst_struc,green,hu,mpi_enreg,paw_dmft,pawang&
737 : !& ,pawprtvol,self_old%qmc_xmu,weiss)
738 :
739 35 : else if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
740 :
741 0 : call ctqmc_calltriqs_c(paw_dmft,green,self_old,hu(:),weiss,self_new,pawprtvol)
742 :
743 35 : else if (abs(paw_dmft%dmft_solv) >= 5) then
744 :
745 : ! == Nothing
746 : ! -------------------
747 : ! call copy_green(weiss,green,opt_tw=1)
748 : ! call copy_green(weiss,green,opt_tw=2)
749 :
750 34 : call qmc_prep_ctqmc(cryst_struc,green,self_old,hu(:),paw_dmft,pawang,pawprtvol,weiss)
751 :
752 :
753 1 : else if (abs(paw_dmft%dmft_solv) == 0) then
754 :
755 : ! == Nothing
756 : ! -------------------
757 : ! weiss%occup%has_operks=0 -> only local part is duplicated
758 1 : call copy_green(weiss,green,opt_tw=2)
759 : end if ! dmft_solv
760 : !call print_green("invWeiss",cryst_struc,weiss,3,paw_dmft,pawtab,2)
761 :
762 : !=======================================================================
763 : !== Treat data from HF QMC
764 : !=======================================================================
765 106 : if (abs(paw_dmft%dmft_solv) >= 4) then
766 : ! propagate qmc_shift (useful for compute_energy)
767 : !if(abs(paw_dmft%dmft_solv)==4) then
768 : ! self_new%qmc_shift(:)=self_old%qmc_shift(:)
769 : ! self_new%qmc_xmu(:)=self_old%qmc_xmu(:)
770 : !end if
771 :
772 : ! == Print local occupations from G(tau)
773 : ! ---------------------------------------
774 :
775 : ! == Fourier back transform of green function G(tau)->G(iw_n) and
776 : ! == compute occupations from g(tau)
777 : ! -------------------------------------------------------------------
778 : !if(abs(paw_dmft%dmft_solv)==4) then
779 : ! write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier Transform t->w of Green Function'
780 : ! call wrtout(std_out,message,'COLL')
781 : ! call fourier_green(cryst_struc,green,paw_dmft,&
782 : !& pawang,opt_ksloc=2,opt_tw=1)
783 : ! do ifreq=1,green%nw
784 : ! xx= green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
785 : ! write(112,*) paw_dmft%omega_lo(ifreq),real(one/xx),aimag(one/xx)
786 : ! write(113,*) paw_dmft%omega_lo(ifreq),real(xx),aimag(xx)
787 : ! end do
788 : ! call flush_unit(112)
789 : ! call flush_unit(113)
790 : ! if(paw_dmft%dmft_solv==5) stop
791 : ! if(pawprtvol>=3) then
792 : ! write(message,'(a,2x,a,f13.5)') ch10,& ! debug
793 : !& " == Print green function for small freq after fourier " ! debug
794 : ! call wrtout(std_out,message,'COLL') ! debug
795 : ! call print_matlu(green%oper(1)%matlu,paw_dmft%natom,1) ! debug
796 : ! end if
797 :
798 : ! write(message,'(2a,i3,13x,a)') ch10,' INVERSE FOURIER OF G0 SUPPRESSED'
799 : ! call wrtout(std_out,message,'COLL')
800 : !end if
801 34 : if (abs(paw_dmft%dmft_solv) == 888) then
802 : ! == Back fourier transform of G_0(tau) for compensation (try or comment or improve FT).
803 : ! -------------------------------------------------------------------
804 0 : write(message,'(2a)') ch10,' === Direct Fourier transform t->w of Weiss'
805 0 : call wrtout(std_out,message,'COLL')
806 0 : call fourier_green(cryst_struc,weiss,paw_dmft,opt_ksloc=2,opt_tw=1)
807 :
808 0 : if (pawprtvol >= 3) then
809 0 : write(message,'(a,2x,a,f13.5)') ch10,& ! debug
810 0 : & " == Print weiss function for small freq after fourier " ! debug
811 0 : call wrtout(std_out,message,'COLL') ! debug
812 0 : call print_matlu(weiss%oper(1)%matlu(:),paw_dmft%natom,1) ! debug
813 : end if ! pawprtvol>=3
814 0 : call destroy_green_tau(weiss)
815 : end if ! dmft_solv=888
816 :
817 : ! == Destroy tau part of green
818 : ! -------------------------------------------------------------------
819 34 : call trace_oper(green%occup_tau,green%charge_ks,green%charge_matlu_solver(:,:),2)
820 34 : green%has_charge_matlu_solver = 2
821 :
822 : end if ! dmft_solv>=5
823 : !=======================================================================
824 : !== End Treat data for QMC
825 : !=======================================================================
826 :
827 : !=======================================================================
828 : !== Integrate green function and printout occupations
829 : !=======================================================================
830 : !For dmft_solv=-1,0,or 1, the green function was not yet computed: it
831 : !cannot be integrated
832 : !=======================================================================
833 106 : if (paw_dmft%dmft_solv >= 2 .and. green%w_type == "imag") then
834 : ! == Integrate G(iw_n)
835 : ! ---------------------
836 103 : write(message,'(2a)') ch10," === Integrate local part of Green's function"
837 103 : call wrtout(std_out,message,'COLL')
838 103 : call integrate_green(green,paw_dmft,2,opt_ksloc=2,opt_after_solver=1)
839 :
840 : ! == Print local occupations from integration of G(iw_n)
841 : ! --------------------------------------------------------
842 103 : call printocc_green(green,5,paw_dmft,3)
843 :
844 : ! == Print G_loc(w)
845 : ! --------------------------------------------------------
846 103 : if (paw_dmft%dmft_prgn == 1) then
847 0 : call print_green('DMFT_IMPURITY',green,1,paw_dmft,opt_wt=1)
848 : end if
849 : end if ! dmft_solv>=2 and w_type="imag"
850 : !stop
851 :
852 : !if(abs(pawprtvol)>0) then
853 : !end if
854 :
855 :
856 : ABI_NVTX_END_RANGE()
857 106 : call timab(622,2,tsec(:))
858 :
859 106 : end subroutine impurity_solve
860 : !!***
861 :
862 : !!****f* ABINIT/dyson
863 : !! NAME
864 : !! dyson
865 : !!
866 : !! FUNCTION
867 : !! Use the Dyson Equation to compute self-energy from green function
868 : !!
869 : !! INPUTS
870 : !! green <type(green_type)>= green function data
871 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
872 : !! self <type(self_type)>= variables related to self-energy
873 : !! weiss <type(green_type)>= Weiss field
874 : !! opt_weissself = 1: compute weiss from green and self
875 : !! = 2: compute self from green and weiss
876 : !!
877 : !! OUTPUT
878 : !!
879 : !! NOTES
880 : !!
881 : !! SOURCE
882 :
883 210 : subroutine dyson(green,paw_dmft,self,weiss,opt_weissself)
884 :
885 : !Arguments ------------------------------------
886 : !scalars
887 : type(green_type), intent(in) :: green
888 : type(paw_dmft_type), intent(in) :: paw_dmft
889 : type(self_type), intent(inout) :: self
890 : type(green_type), intent(inout) :: weiss
891 : integer, intent(in) :: opt_weissself
892 : ! type(paw_dmft_type), intent(inout) :: paw_dmft
893 : !Local variables ------------------------------
894 : integer :: ifreq,myproc,natom,nspinor,nsppol,weissinv
895 : logical :: triqs
896 : real(dp) :: tsec(2)
897 210 : type(matlu_type), allocatable :: greeninv(:)
898 : character(len=500) :: message
899 : ! type
900 : ! type(matlu_type), pointer :: matlutemp,matlu1,matlu2
901 : !************************************************************************
902 :
903 210 : call timab(623,1,tsec(:))
904 : DBG_ENTER("COLL")
905 :
906 210 : myproc = paw_dmft%myproc
907 210 : natom = paw_dmft%natom
908 210 : nsppol = paw_dmft%nsppol
909 210 : nspinor = paw_dmft%nspinor
910 210 : triqs = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
911 210 : weissinv = merge(0,1,paw_dmft%dmft_solv==2.or.triqs)
912 :
913 210 : if (opt_weissself == 1) then
914 106 : write(message,'(2a)') ch10," === Use Dyson's Equation => weiss"
915 106 : call wrtout(std_out,message,'COLL')
916 104 : else if (opt_weissself == 2) then
917 104 : write(message,'(2a)') ch10," === Use Dyson's Equation => self"
918 104 : call wrtout(std_out,message,'COLL')
919 : end if ! opt_weisself
920 :
921 : !call xmpi_barrier(spaceComm)
922 :
923 1366 : ABI_MALLOC(greeninv,(natom))
924 210 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),greeninv(:))
925 :
926 46324 : do ifreq=1,green%nw
927 :
928 46114 : if (green%distrib%procf(ifreq) /= myproc) cycle
929 :
930 16000 : call copy_matlu(green%oper(ifreq)%matlu(:),greeninv(:),natom)
931 16000 : call inverse_matlu(greeninv(:),natom)
932 :
933 16210 : if (opt_weissself == 1) then
934 :
935 : ! warning green is now inversed
936 8032 : call add_matlu(greeninv(:),self%oper(ifreq)%matlu(:),weiss%oper(ifreq)%matlu(:),natom,1)
937 8032 : if (.not. triqs) then
938 8032 : call inverse_oper(weiss%oper(ifreq),2)
939 : end if
940 :
941 7968 : else if (opt_weissself == 2) then
942 :
943 : ! write(59,*) paw_dmft%omega_lo(ifreq), real(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
944 : ! write(61,*) paw_dmft%omega_lo(ifreq), real(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
945 : ! write(60,*) paw_dmft%omega_lo(ifreq), real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
946 : ! call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1)
947 : ! write(62,*) paw_dmft%omega_lo(ifreq), real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
948 : ! call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1)
949 : ! write(63,*) paw_dmft%omega_lo(ifreq), real(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
950 :
951 : ! write(std_out,*) "-----------------------IFREQ",ifreq
952 : ! call print_matlu(greeninv%oper(ifreq)%matlu,paw_dmft%natom,1,opt_diag=-1)
953 : !call inverse_oper(greeninv%oper(ifreq),option=1,prtopt=1)
954 : ! call print_matlu(greeninv%oper(ifreq)%matlu,paw_dmft%natom,1,opt_diag=-1)
955 : ! If paw_dmft%dmft_solv==2, then inverse of weiss function is
956 : ! computed in m_hubbard_one.F90
957 7968 : if (weissinv /= 0) then
958 6288 : call inverse_oper(weiss%oper(ifreq),2)
959 : end if
960 :
961 : ! write(std_out,*) weiss%oper(1)%matlu(ifreq)%mat(1,1,1,1,1),"-",greeninv%oper(ifreq)
962 7968 : call add_matlu(weiss%oper(ifreq)%matlu(:),greeninv(:),self%oper(ifreq)%matlu(:),natom,-1)
963 :
964 : ! write(64,*) paw_dmft%omega_lo(ifreq), real(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
965 : ! write(65,*) paw_dmft%omega_lo(ifreq), real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
966 : ! write(66,*) paw_dmft%omega_lo(ifreq), real(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
967 : else
968 0 : message = " BUG in dyson.F90"
969 0 : ABI_BUG(message)
970 : end if ! opt_weissself
971 : end do ! ifreq
972 :
973 210 : if (opt_weissself == 1) then
974 106 : call gather_oper(weiss%oper(:),weiss%distrib,paw_dmft,opt_ksloc=2)
975 104 : else if (opt_weissself == 2) then
976 104 : call gather_oper(self%oper(:),self%distrib,paw_dmft,opt_ksloc=2)
977 : end if ! opt_weisself
978 :
979 210 : call destroy_matlu(greeninv(:),natom)
980 946 : ABI_FREE(greeninv)
981 :
982 210 : call timab(623,2,tsec(:))
983 : DBG_EXIT("COLL")
984 :
985 210 : end subroutine dyson
986 : !!***
987 :
988 : !!****f* m_dmft/spectral_function
989 : !! NAME
990 : !! spectral_function
991 : !!
992 : !! FUNCTION
993 : !! Print the spectral function computed from Green function in real frequency
994 : !!
995 : !! INPUTS
996 : !! cryst_struc <type(crystal_t)>=crystal structure data
997 : !! green <type(green_type)>= green function data
998 : !! hu <type(hu_type)>= datatype of type hu
999 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
1000 : !! self <type(self_type)>= variables related to self-energy
1001 : !! prtopt= option for printing
1002 : !!
1003 : !! OUTPUT
1004 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
1005 : !!
1006 : !! NOTES
1007 : !!
1008 : !! SOURCE
1009 :
1010 0 : subroutine spectral_function(cryst_struc,green,hu,paw_dmft,&
1011 0 : & pawtab,self_old,prtopt)
1012 :
1013 : use m_dftu_self, only : dftu_self
1014 : use m_green, only : compute_green,copy_green,destroy_green,init_green,print_green
1015 : use m_hubbard_one, only : hubbard_one
1016 : use m_matlu, only : copy_matlu
1017 : use m_self, only : dc_self,destroy_self,initialize_self,rw_self
1018 :
1019 : !Arguments ------------------------------------
1020 : !scalars
1021 : type(crystal_t),intent(in) :: cryst_struc
1022 : type(green_type), intent(in) :: green
1023 : type(hu_type),intent(inout) :: hu(cryst_struc%ntypat)
1024 : !type(MPI_type), intent(inout) :: mpi_enreg
1025 : type(pawtab_type),intent(inout) :: pawtab(cryst_struc%ntypat)
1026 : type(self_type), intent(inout) :: self_old
1027 : type(paw_dmft_type), intent(inout) :: paw_dmft
1028 : integer, intent(in) :: prtopt
1029 :
1030 : !Local variables ------------------------------
1031 : character(len=500) :: message
1032 0 : type(green_type) :: greenr
1033 0 : type(green_type) :: weissr
1034 0 : type(self_type) :: selfr
1035 : !scalars
1036 : !************************************************************************
1037 : !character(len=500) :: message
1038 :
1039 : ! opt_oper_ksloc=3 to be able to compute spectral function.
1040 0 : call init_green(greenr,paw_dmft,opt_oper_ksloc=3,wtype="real")
1041 0 : call init_green(weissr,paw_dmft,wtype="real")
1042 0 : call copy_matlu(green%occup%matlu,greenr%occup%matlu,paw_dmft%natom)
1043 0 : call initialize_self(selfr,paw_dmft,wtype="real")
1044 : !=======================================================================
1045 : !== Solve impurity model with green function for real frequency
1046 : !=======================================================================
1047 0 : write(message,'(2a,i3,13x,a)') ch10,' === Write Spectral function'
1048 0 : call wrtout(std_out,message,'COLL')
1049 0 : if(abs(paw_dmft%dmft_solv)==1) then
1050 :
1051 : ! == DFT+U for test
1052 : ! -------------------
1053 : call dftu_self(cryst_struc,greenr,paw_dmft,&
1054 0 : & pawtab,selfr,opt_dftu=1,prtopt=prtopt)
1055 0 : else if(abs(paw_dmft%dmft_solv)==2) then
1056 :
1057 : ! == Hubbard One
1058 : ! -------------------
1059 : call hubbard_one(cryst_struc,greenr,hu,paw_dmft,&
1060 0 : & prtopt,self_old%hdc,weissr)
1061 :
1062 0 : else if(abs(paw_dmft%dmft_solv)==4) then
1063 :
1064 : ! == Nothing
1065 : ! -------------------
1066 0 : message = "spectral_function: This section of code is disabled!"
1067 0 : ABI_ERROR(message)
1068 0 : call copy_green(weissr,greenr,opt_tw=1)
1069 :
1070 0 : else if(abs(paw_dmft%dmft_solv)>=5) then
1071 :
1072 : ! == Nothing
1073 : ! -------------------
1074 0 : ABI_ERROR("Stopping before copy_green")
1075 0 : call copy_green(weissr,greenr,opt_tw=1)
1076 :
1077 0 : else if(abs(paw_dmft%dmft_solv)==0) then
1078 :
1079 : ! == Nothing
1080 : ! -------------------
1081 : ! weiss%occup%has_operks=0 -> only local part is duplicated
1082 0 : call copy_green(weissr,greenr,opt_tw=2)
1083 : end if
1084 :
1085 : !=======================================================================
1086 : !== Integrate green function and printout occupations
1087 : !For dmft_solv=-1,0,or 1 , the green function was not computed: it
1088 : !cannot be integrated
1089 : !=======================================================================
1090 0 : call dc_self(green%charge_matlu_solver,selfr%hdc%matlu,hu,paw_dmft,pawtab,green%occup%matlu)
1091 0 : if(abs(paw_dmft%dmft_solv)/=1.and.paw_dmft%dmft_solv/=0) then
1092 0 : call dyson(greenr,paw_dmft,selfr,weissr,opt_weissself=2)
1093 : end if
1094 0 : call compute_green(greenr,paw_dmft,1,selfr,opt_self=1)
1095 0 : call print_green("realw",greenr,4,paw_dmft)
1096 0 : call rw_self(selfr,paw_dmft,prtopt=2,opt_rw=2)
1097 :
1098 : if(abs(prtopt)>0) then
1099 : end if
1100 0 : call destroy_self(selfr)
1101 0 : call destroy_green(weissr)
1102 0 : call destroy_green(greenr)
1103 :
1104 0 : end subroutine spectral_function
1105 : !!***
1106 :
1107 : END MODULE m_dmft
1108 : !!***
|
