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