Line data Source code
1 : !!****m* ABINIT/m_energy
2 : !! NAME
3 : !! m_energy
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 :
24 : #include "abi_common.h"
25 :
26 : MODULE m_energy
27 :
28 : use defs_basis
29 : use m_abi_linalg, only : abi_xgemm
30 : use m_abicore
31 : use m_errors
32 : use m_green, only : green_type,occup_fd
33 : use m_matlu, only : add_matlu,copy_matlu,destroy_matlu,init_matlu, &
34 : & matlu_type,trace_prod_matlu
35 : use m_oper, only : oper_type
36 : use m_paw_correlations, only : pawuenergy
37 : use m_paw_dmft, only : paw_dmft_type
38 : use m_pawtab, only : pawtab_type
39 : use m_self, only : self_type
40 : use m_xmpi, only : xmpi_sum
41 :
42 : implicit none
43 :
44 : private
45 :
46 : public :: init_energy
47 : public :: compute_energy
48 : public :: compute_migdal_energy
49 : public :: compute_dftu_energy
50 : public :: destroy_energy
51 : public :: print_energy
52 : public :: compute_noninterentropy
53 : public :: compute_free_energy
54 : public :: compute_trace_log_loc
55 : public :: print_free_energy
56 : !!***
57 :
58 : !!****t* m_energy/energy_type
59 : !! NAME
60 : !! energy_type
61 : !!
62 : !! FUNCTION
63 : !! This structured datatype contains interaction matrices for the correlated subspace
64 : !!
65 : !! SOURCE
66 :
67 : type, public :: energy_type ! for each typat
68 :
69 : real(dp) :: e_dc_tot
70 :
71 : real(dp) :: e_dcdc
72 :
73 : real(dp) :: e_hu_dftu_tot
74 :
75 : real(dp) :: e_hu_mig_tot
76 :
77 : real(dp) :: e_hu_qmc_tot
78 :
79 : real(dp) :: e_hu_tot
80 :
81 : real(dp) :: eband_dft
82 :
83 : real(dp) :: eband_dmft
84 :
85 : real(dp) :: edmft
86 :
87 : real(dp) :: ekin_imp
88 :
89 : real(dp) :: emig_imp
90 :
91 : real(dp) :: emig_loc
92 :
93 : real(dp) :: fband_dft
94 :
95 : real(dp) :: fband_dmft
96 :
97 : real(dp) :: fband_imp
98 :
99 : real(dp) :: fband_weiss
100 :
101 : real(dp) :: fdmft
102 :
103 : real(dp) :: fimp
104 :
105 : real(dp) :: integral
106 :
107 : !real(dp) :: natom
108 :
109 : real(dp) :: sdmft
110 :
111 : real(dp) :: simp
112 :
113 : real(dp), allocatable :: e_dc(:)
114 :
115 : real(dp), allocatable :: e_hu_dftu(:)
116 :
117 : real(dp), allocatable :: e_hu_mig(:)
118 :
119 : real(dp), allocatable :: e_hu_qmc(:)
120 :
121 : real(dp), ABI_CONTIGUOUS pointer :: e_hu(:) => null()
122 :
123 : end type energy_type
124 :
125 : !!***
126 :
127 : !----------------------------------------------------------------------
128 :
129 :
130 : CONTAINS !========================================================================================
131 : !!***
132 :
133 : !!****f* m_energy/init_energy
134 : !! NAME
135 : !! init_energy
136 : !!
137 : !! FUNCTION
138 : !! Allocate variables used in type energy_type.
139 : !!
140 : !! INPUTS
141 : !!
142 : !! OUTPUTS
143 : !! energies_dmft = datastructure for dmft energy
144 : !! natom = number of atoms
145 : !!
146 : !! SOURCE
147 :
148 102 : subroutine init_energy(energies_dmft,natom)
149 :
150 : !Arguments ------------------------------------
151 : !type
152 : type(energy_type), target, intent(inout) :: energies_dmft
153 : integer, intent(in) :: natom
154 : !Local variables ------------------------------------
155 : !************************************************************************
156 :
157 306 : ABI_MALLOC(energies_dmft%e_dc,(natom))
158 204 : ABI_MALLOC(energies_dmft%e_hu_dftu,(natom))
159 204 : ABI_MALLOC(energies_dmft%e_hu_mig,(natom))
160 204 : ABI_MALLOC(energies_dmft%e_hu_qmc,(natom))
161 102 : energies_dmft%e_hu => energies_dmft%e_hu_mig(:)
162 456 : energies_dmft%e_dc(:) = zero
163 456 : energies_dmft%e_hu_dftu(:) = zero
164 456 : energies_dmft%e_hu_mig(:) = zero
165 456 : energies_dmft%e_hu_qmc(:) = zero
166 102 : energies_dmft%e_dc_tot = zero
167 102 : energies_dmft%e_dcdc = zero
168 102 : energies_dmft%e_hu_dftu_tot = zero
169 102 : energies_dmft%e_hu_mig_tot = zero
170 102 : energies_dmft%e_hu_qmc_tot = zero
171 102 : energies_dmft%e_hu_tot = zero
172 102 : energies_dmft%eband_dft = zero
173 102 : energies_dmft%eband_dmft = zero
174 102 : energies_dmft%edmft = zero
175 102 : energies_dmft%ekin_imp = zero
176 102 : energies_dmft%emig_imp = zero
177 102 : energies_dmft%emig_loc = zero
178 102 : energies_dmft%fband_dft = zero
179 102 : energies_dmft%fband_dmft = zero
180 102 : energies_dmft%fband_imp = zero
181 102 : energies_dmft%fband_weiss = zero
182 102 : energies_dmft%fdmft = zero
183 102 : energies_dmft%fimp = zero
184 102 : energies_dmft%integral = zero
185 : !energies_dmft%natom = natom
186 102 : energies_dmft%sdmft = zero
187 102 : energies_dmft%simp = zero
188 :
189 102 : end subroutine init_energy
190 : !!***
191 :
192 : !!****f* m_energy/destroy_energy
193 : !! NAME
194 : !! destroy_energy
195 : !!
196 : !! FUNCTION
197 : !! Deallocate energies_dmft
198 : !!
199 : !! INPUTS
200 : !! energies_dmft = datastructure for dmft energy
201 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
202 : !!
203 : !! OUTPUT
204 : !!
205 : !! SOURCE
206 :
207 102 : subroutine destroy_energy(energies_dmft,paw_dmft)
208 :
209 : !Arguments ------------------------------------
210 : !scalars
211 : type(energy_type), intent(inout) :: energies_dmft
212 : type(paw_dmft_type), intent(inout) :: paw_dmft
213 : !Local variables-------------------------------
214 : ! *********************************************************************
215 :
216 102 : paw_dmft%e_dc = energies_dmft%e_dc_tot
217 102 : paw_dmft%e_hu = energies_dmft%e_hu_tot
218 102 : paw_dmft%sdmft = energies_dmft%sdmft
219 102 : paw_dmft%simp = energies_dmft%simp
220 :
221 102 : energies_dmft%e_hu => null()
222 102 : ABI_SFREE(energies_dmft%e_dc)
223 102 : ABI_SFREE(energies_dmft%e_hu_dftu)
224 102 : ABI_SFREE(energies_dmft%e_hu_mig)
225 102 : ABI_SFREE(energies_dmft%e_hu_qmc)
226 :
227 102 : end subroutine destroy_energy
228 : !!***
229 :
230 : !!****f* m_energy/print_energy
231 : !! NAME
232 : !! print_energy
233 : !!
234 : !! FUNCTION
235 : !! Print different components of DMFT contribution to the internal energy.
236 : !!
237 : !! INPUTS
238 : !! energies_dmft = datastructure for dmft energy
239 : !! pawprtvol = flag for print_energy
240 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
241 : !! idmftloop = iteration number of the DFT+DMFT loop
242 : !!
243 : !! OUTPUT
244 : !!
245 : !! SIDE EFFECTS
246 : !!
247 : !! SOURCE
248 :
249 292 : subroutine print_energy(energies_dmft,pawprtvol,paw_dmft,idmftloop)
250 :
251 : !Arguments ------------------------------------
252 : !type
253 : type(energy_type), intent(in) :: energies_dmft
254 : type(paw_dmft_type), intent(in) :: paw_dmft
255 : integer, intent(in) :: pawprtvol,idmftloop
256 : !Local variables-------------------------------
257 : integer :: iatom,lpawu
258 : character(len=4) :: tag
259 : character(len=1000) :: message
260 : ! *********************************************************************
261 :
262 292 : if (abs(pawprtvol) >= 3) then
263 908 : do iatom=1,paw_dmft%natom
264 700 : lpawu = paw_dmft%lpawu(iatom)
265 700 : if (lpawu == -1) cycle
266 208 : write(tag,'(i4)') iatom
267 208 : write(message,'(a,4x,3a)') ch10,"For Correlated Atom ",trim(adjustl(tag)),","
268 208 : call wrtout(std_out,message,'COLL')
269 208 : write(message,'(26x,a,f12.6)') " E_hu =",energies_dmft%e_hu(iatom)
270 208 : call wrtout(std_out,message,'COLL')
271 208 : write(message,'(26x,a,f12.6)') " E_hu_mig =",energies_dmft%e_hu_mig(iatom)
272 208 : call wrtout(std_out,message,'COLL')
273 208 : write(message,'(26x,a,f12.6)') " E_hu_qmc =",energies_dmft%e_hu_qmc(iatom)
274 208 : call wrtout(std_out,message,'COLL')
275 208 : write(message,'(26x,a,f12.6)') " E_hu_dftu =",energies_dmft%e_hu_dftu(iatom)
276 208 : call wrtout(std_out,message,'COLL')
277 208 : write(message,'(26x,a,f12.6)') " E_dc =",energies_dmft%e_dc(iatom)
278 908 : call wrtout(std_out,message,'COLL')
279 : end do ! iatom
280 : end if ! abs(pawprtvol)>=3
281 292 : write(message,'(a,5x,2a,5x,a,9(a,5x,a,2x,f15.11),a,5x,a)') ch10, &
282 292 : & "-----------------------------------------------",ch10, &
283 292 : & "--- Energy in DMFT (in Ha) ",ch10, &
284 292 : & "--- E_bandlda (1) (Ha.) = ",energies_dmft%eband_dft,ch10, &
285 292 : & "--- E_banddmft(2) (Ha.) = ",energies_dmft%eband_dmft,ch10, &
286 292 : & "--- E_hu (3) (Ha.) = ",energies_dmft%e_hu_tot,ch10, &
287 292 : & "--- E_hu_mig (4) (Ha.) = ",energies_dmft%e_hu_mig_tot,ch10, &
288 292 : & "--- E_hu_qmc (4) (Ha.) = ",energies_dmft%e_hu_qmc_tot,ch10, &
289 292 : & "--- E_hu_dftu (5) (Ha.) = ",energies_dmft%e_hu_dftu_tot,ch10, &
290 292 : & "--- E_dc (6) (Ha.) = ",energies_dmft%e_dc_tot,ch10, &
291 292 : & "--- edmft=( 3-6)(Ha.) = ",energies_dmft%edmft,ch10, &
292 292 : & "--- (2-1+3-6)(Ha.) = ",energies_dmft%eband_dmft-energies_dmft%eband_dft+energies_dmft%edmft,ch10, &
293 584 : & "-----------------------------------------------"
294 292 : call wrtout(std_out,message,'COLL')
295 292 : if (idmftloop >= 1) then
296 106 : write(message,'(a,i3,1x,f15.11,a)') " (Edmft",idmftloop,energies_dmft%edmft,")"
297 106 : call wrtout(ab_out,message,'COLL')
298 : end if
299 :
300 292 : end subroutine print_energy
301 : !!***
302 :
303 : !!****f* m_energy/compute_energy
304 : !! NAME
305 : !! compute_energy
306 : !!
307 : !! FUNCTION
308 : !! Compute and print the different contributions for the DMFT energy.
309 : !!
310 : !! INPUTS
311 : !! energies_dmft = datastructure for dmft energy
312 : !! green <type(green_type)>= green function data
313 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
314 : !! pawprtvol = flag for printing
315 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
316 : !! self <type(self_type)>= self energy function data
317 : !! occ_type= character ("lda" or "nlda") for printing.
318 : !! part = "band" : compute the DFT and DMFT band energies
319 : !! = "corr" : compute the interaction energy
320 : !! = "both" : compute band energy and interaction energy
321 : !! = "none" : do not print
322 : !!
323 : !! OUTPUT
324 : !!
325 : !! SIDE EFFECTS
326 : !! energies_dmft = datastructure for dmft energy
327 : !!
328 : !! SOURCE
329 :
330 398 : subroutine compute_energy(energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type,part)
331 :
332 : !Arguments ------------------------------------
333 : !type
334 : type(energy_type), target, intent(inout) :: energies_dmft
335 : type(green_type), intent(in) :: green
336 : type(paw_dmft_type), intent(in) :: paw_dmft
337 : type(pawtab_type), intent(in) :: pawtab(paw_dmft%ntypat)
338 : type(self_type), intent(in) :: self
339 : integer, intent(in) :: pawprtvol
340 : character(len=4), intent(in) :: occ_type,part
341 : ! integer :: prtopt
342 : !Local variables-------------------------------
343 : integer :: iatom,lpawu
344 : character(len=500) :: message
345 : ! *********************************************************************
346 :
347 398 : if (part == 'both') then
348 199 : if (occ_type == " lda") then
349 93 : write(message,'(2a)') ch10," == Check: Compute DFT energy terms"
350 : else
351 106 : write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms"
352 : end if
353 199 : call wrtout(std_out,message,'COLL')
354 199 : else if (part == 'band') then
355 93 : write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms : Band energy terms"
356 93 : call wrtout(std_out,message,'COLL')
357 106 : else if (part == 'corr') then
358 0 : write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms : Correlation energy terms only"
359 0 : call wrtout(std_out,message,'COLL')
360 : !else if(part=='none') then
361 : end if ! part
362 :
363 : ! Only imaginary frequencies here
364 398 : if (green%w_type == "real" .or. self%w_type == "real") then
365 0 : message = 'compute_energy not implemented for real frequency'
366 0 : ABI_BUG(message)
367 : end if
368 :
369 398 : if (part == 'band' .or. part == 'both') then
370 292 : call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft=1)
371 : end if
372 : ! == Compute Band Energy Alternative version: two steps
373 : ! == Compute Tr[ln G^{-1}] and -Tr[(Self-hdc)G_dmft]
374 : ! -----------------------------------------------------------------------
375 : ! if (part == 'band') then ! ie if thdyn="fcalc" in m_dmft.F90
376 : ! call compute_B3(cryst_struc,energies_dmft,eband2,green,mpi_enreg,paw_dmft,2,pawang,self,occ_type,0)
377 : ! write(message,'(2a,f10.6)') ch10,"Compute Band energy test KS ",eband2
378 : ! call wrtout(std_out,message,'COLL')
379 : ! call wrtout(ab_out,message,'COLL')
380 : !
381 : ! call compute_B3(cryst_struc,energies_dmft,eband2,green,mpi_enreg,paw_dmft,2,pawang,self,occ_type,1)
382 : ! write(message,'(2a,f10.6)') ch10,"Compute Band energy test Self statique",eband2
383 : ! call wrtout(std_out,message,'COLL')
384 : ! call wrtout(ab_out,message,'COLL')
385 : !
386 : !! == Compute Band Energy (classical)
387 : !! -----------------------------------------------------------------------
388 : ! call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,fcalc_dft=3)
389 : ! write(std_out,*) paw_dmft%fermie_dft,paw_dmft%fermie
390 : ! write(message,'(2a,f10.6)') ch10,"Compute Band energy ref free lda -ef ",energies_dmft%eband_dft
391 : ! call wrtout(std_out,message,'COLL')
392 : ! call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft=1)
393 : ! write(message,'(2a,f10.6)') ch10,"Compute Band energy ref -ef ",energies_dmft%eband_dmft
394 : ! call wrtout(std_out,message,'COLL')
395 : !! call wrtout(ab_out,message,'COLL')
396 : ! write(message,'(2a,f10.6)') ch10,"Compute Band energy ref lda ",energies_dmft%eband_dft
397 : ! call wrtout(std_out,message,'COLL')
398 : !! if(occ_type=="nlda") eband2=energies_dmft%eband_dmft
399 : ! else
400 : !call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft=1)
401 : !endif
402 :
403 : !end if
404 :
405 398 : if (part == 'corr' .or. part == 'both') then
406 :
407 : ! == Compute Correlation energy from Migdal formula
408 : ! -----------------------------------------------------------------------
409 199 : if (occ_type /= " lda") then
410 106 : call compute_migdal_energy(energies_dmft%e_hu_mig(:),energies_dmft%e_hu_mig_tot,green,paw_dmft,self)
411 : end if
412 : ! write(std_out,*) "MIGDAL",e_hu_migdal_tot,e_hu_migdal
413 :
414 : ! == Compute DFT+U interaction energy
415 : ! -----------------------------------------------------------------------
416 199 : call compute_dftu_energy(energies_dmft,green,paw_dmft,pawtab(:))
417 199 : if (abs(paw_dmft%dmft_solv) <= 1) then
418 6 : energies_dmft%e_hu => energies_dmft%e_hu_dftu(:)
419 6 : energies_dmft%e_hu_tot = energies_dmft%e_hu_dftu_tot
420 6 : if ((abs(energies_dmft%e_hu_tot-energies_dmft%e_hu_mig_tot) >= tol6) .and. (occ_type /= " lda")) then
421 0 : write(message,'(2a,2e18.8,2x,a)') ch10,' BUG: Migdal energy and DFT+U energy do not coincide',&
422 0 : & energies_dmft%e_hu_tot,energies_dmft%e_hu_mig_tot,occ_type
423 0 : ABI_ERROR(message)
424 : end if
425 : else if (paw_dmft%dmft_solv == 2 .or. ((paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) &
426 193 : & .and. (.not. paw_dmft%dmft_triqs_measure_density_matrix)) .or. paw_dmft%dmft_solv == 9) then
427 125 : energies_dmft%e_hu => energies_dmft%e_hu_mig(:)
428 125 : energies_dmft%e_hu_tot = energies_dmft%e_hu_mig_tot
429 125 : energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_tot
430 68 : else if (paw_dmft%dmft_solv == 5 .or. paw_dmft%dmft_solv == 8 .or. paw_dmft%dmft_solv == 10 .or.&
431 : & ((paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) .and. &
432 : & paw_dmft%dmft_triqs_measure_density_matrix) .and. occ_type /= " lda") then
433 68 : if (paw_dmft%dmft_solv == 8 .or. paw_dmft%dmft_solv == 10) then
434 0 : write(message,'(2a)') ch10,"Warning, energy is recently computed, not checked"
435 0 : call wrtout(std_out,message,'COLL')
436 : end if
437 : ! == Compute Correlation energy from QMC correlations.
438 : ! -----------------------------------------------------------------------
439 68 : energies_dmft%e_hu_qmc_tot = zero
440 192 : do iatom=1,paw_dmft%natom
441 124 : lpawu = paw_dmft%lpawu(iatom)
442 124 : if (lpawu == -1) cycle
443 68 : energies_dmft%e_hu_qmc(iatom) = green%ecorr_qmc(iatom)
444 192 : energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_qmc_tot + energies_dmft%e_hu_qmc(iatom)
445 : end do ! iatom
446 68 : energies_dmft%e_hu => energies_dmft%e_hu_qmc(:)
447 68 : energies_dmft%e_hu_tot = energies_dmft%e_hu_qmc_tot
448 : end if ! dmft_solv
449 : ! energies_dmft%edmft=energies_dmft%e_hu_mig_tot-energies_dmft%e_dc_tot
450 199 : energies_dmft%edmft = energies_dmft%e_hu_tot - energies_dmft%e_dc_tot
451 :
452 : end if ! part
453 :
454 : ! if(part='corr'.or.part='both') then
455 398 : if (part /= 'none') then
456 292 : call print_energy(energies_dmft,pawprtvol,paw_dmft,paw_dmft%idmftloop)
457 : end if
458 : ! write(message,'(2a)') ch10," == The DFT+U self-energy is == "
459 : ! call wrtout(std_out,message,'COLL')
460 : ! call print_oper(self%oper(1),5,paw_dmft,2)
461 : ! a voir: energies_dmft%e_hu_tot = energies_dmft%e_hu_dftu_tot
462 :
463 398 : end subroutine compute_energy
464 : !!***
465 :
466 : !!****f* m_energy/compute_band_energy
467 : !! NAME
468 : !! compute_band_energy
469 : !!
470 : !! FUNCTION
471 : !! Compute the DFT and DMFT band energy
472 : !!
473 : !! INPUTS
474 : !! energies_dmft = datastructure for dmft energy
475 : !! green <type(green_type)>= green function data
476 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
477 : !! occ_type= character ("lda" or "nlda") for printing.
478 : !! fcalc_dft= if present, compute free energy/grand potential instead of total energy.
479 : !! = 2/4 use the DMFT Fermi level
480 : !! = 2/3 compute grand potential
481 : !! = 1/3 use the DFT Fermi level
482 : !! ecalc_dft= 2/4 use the DMFT Fermi level
483 : !! = 2/3 substract mu*N to the DFT band energy
484 : !! = 1/3 use the DFT Fermi level
485 : !! ecalc_dmft= if present, substract mu*N to the DMFT band energy
486 : !!
487 : !! OUTPUT
488 : !!
489 : !! SIDE EFFECTS
490 : !! energies_dmft = datastructure for dmft energy
491 : !!
492 : !! SOURCE
493 :
494 292 : subroutine compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft,fcalc_dft,ecalc_dmft)
495 :
496 : !Arguments ------------------------------------
497 : type(energy_type), intent(inout) :: energies_dmft
498 : type(green_type), intent(in) :: green
499 : type(paw_dmft_type), intent(in) :: paw_dmft
500 : character(len=4), intent(in) :: occ_type
501 : integer, optional, intent(in) :: ecalc_dft,ecalc_dmft,fcalc_dft
502 : ! integer :: prtopt
503 : !Local variables-------------------------------
504 : integer :: band_index,ib,ibc,ikpt,isppol,nband_k,nkpt,nspinor,nsppol
505 : real(dp) :: beta,eig,fermie_used,occ,totch2,wtk !,totch3
506 : character(len=500) :: message
507 : ! *********************************************************************
508 :
509 292 : if (occ_type == " lda") then
510 93 : write(message,'(2a)') ch10," == Compute DFT Band Energy terms"
511 199 : else if (present(fcalc_dft)) then
512 0 : write(message,'(2a)') ch10," == Compute DFT Free Energy terms"
513 : else
514 199 : write(message,'(2a)') ch10," == Compute DMFT Band Energy terms"
515 : end if
516 292 : call wrtout(std_out,message,'COLL')
517 292 : beta = one / paw_dmft%temp
518 :
519 : ! == Compute Band Energy
520 : ! -----------------------------------------------------------------------
521 292 : if (occ_type == " lda") energies_dmft%eband_dft = zero
522 292 : if (.not. present(fcalc_dft)) energies_dmft%eband_dmft = zero
523 0 : if (present(fcalc_dft)) energies_dmft%fband_dft = zero
524 : !totch = zero
525 292 : totch2 = zero
526 : !totch3 = zero
527 :
528 292 : nkpt = paw_dmft%nkpt
529 292 : nspinor = paw_dmft%nspinor
530 292 : nsppol = paw_dmft%nsppol
531 :
532 292 : band_index = 0
533 710 : do isppol=1,nsppol
534 4010 : do ikpt=1,nkpt
535 3300 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
536 3300 : wtk = paw_dmft%wtk(ikpt)
537 3300 : ibc = 0
538 80664 : do ib=1,nband_k
539 77364 : if ((.not. paw_dmft%band_in(ib)) .and. (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7)) cycle
540 31268 : if (paw_dmft%band_in(ib)) ibc = ibc + 1
541 31268 : eig = paw_dmft%eigen(ib+band_index)
542 34568 : if (present(fcalc_dft)) then
543 0 : if (fcalc_dft == 1 .or. fcalc_dft == 3) fermie_used = paw_dmft%fermie_dft
544 0 : if (fcalc_dft == 2 .or. fcalc_dft == 4) fermie_used = paw_dmft%fermie ! only for B3 terms
545 : energies_dmft%fband_dft = energies_dmft%fband_dft + wtk*merge(-log(one+exp(-beta*(eig-fermie_used))), &
546 0 : & (beta*(eig-fermie_used)-log(one+exp(beta*(eig-fermie_used)))),(eig-fermie_used)>=zero)
547 : else
548 31268 : if (occ_type == " lda") then ! usual calculation: total non interacting energy
549 9852 : fermie_used = paw_dmft%fermie_dft
550 : ! write(std_out,*) "isppol,ikpt,ib",isppol,ikpt,ib
551 : ! write(std_out,*) "paw_dmft%eigen_dft",paw_dmft%eigen_dft(isppol,ikpt,ib)
552 : ! write(std_out,*) green%occup%ks(isppol,ikpt,ib,ib)
553 : ! write(std_out,*) occup_fd(paw_dmft%eigen_dft(isppol,ikpt,ib),paw_dmft%fermie,paw_dmft%temp)
554 9852 : if (present(ecalc_dft)) then
555 9852 : if (ecalc_dft == 1 .or. ecalc_dft == 3) fermie_used = paw_dmft%fermie_dft
556 9852 : if (ecalc_dft == 2 .or. ecalc_dft == 4) fermie_used = paw_dmft%fermie ! only for B3 terms
557 9852 : occ = occup_fd(eig,fermie_used,paw_dmft%temp)
558 9852 : if (ecalc_dft == 3 .or. ecalc_dft == 2) then
559 0 : energies_dmft%eband_dft = energies_dmft%eband_dft - occ*fermie_used*wtk
560 0 : totch2 = totch2 + wtk*occ
561 : end if
562 : else
563 0 : occ = occup_fd(eig,fermie_used,paw_dmft%temp)
564 : end if ! present(ecalc_dft)
565 9852 : energies_dmft%eband_dft = energies_dmft%eband_dft + occ*eig*wtk
566 : end if ! occ_type=" lda"
567 31268 : if (paw_dmft%band_in(ib)) then
568 31268 : occ = dble(green%occup%ks(ibc,ibc,ikpt,isppol))
569 : else
570 0 : occ = occup_fd(eig,paw_dmft%fermie,paw_dmft%temp)
571 : end if
572 31268 : energies_dmft%eband_dmft = energies_dmft%eband_dmft + occ*eig*wtk
573 : !totch3 = totch3 + paw_dmft%wtk(ikpt)*green%occup%ks(ib,ib,ikpt,isppol)
574 31268 : if (present(ecalc_dmft)) energies_dmft%eband_dmft = energies_dmft%eband_dmft - &
575 0 : & occ*paw_dmft%fermie*wtk
576 : end if ! present(fcalc_dft)
577 : !end if
578 : end do ! ib
579 3718 : band_index = band_index + nband_k
580 : end do ! ikpt
581 : end do ! isppol
582 :
583 292 : if (present(fcalc_dft)) then
584 0 : energies_dmft%fband_dft = energies_dmft%fband_dft * paw_dmft%temp
585 0 : if (nsppol == 1 .and. nspinor == 1) energies_dmft%fband_dft = energies_dmft%fband_dft * two
586 0 : if (fcalc_dft == 1 .or. fcalc_dft == 4) energies_dmft%fband_dft = energies_dmft%fband_dft + &
587 0 : & fermie_used*paw_dmft%nelectval
588 : else
589 292 : if (occ_type == " lda" .and. nsppol == 1 .and. nspinor == 1) energies_dmft%eband_dft = two * energies_dmft%eband_dft
590 292 : if (nsppol == 1 .and. nspinor == 1) energies_dmft%eband_dmft = two * energies_dmft%eband_dmft
591 : end if ! present(fcalc_dft)
592 : !if (fcalc_dft == 3 .or. fcalc_dft == 2) write(std_out,*) "compute_band_energy totch",totch
593 : !end if
594 :
595 292 : if (present(ecalc_dft)) then
596 292 : if (ecalc_dft == 3 .or. ecalc_dft == 2) write(std_out,*) "compute_band_energy totch2",totch2
597 : end if
598 : ! write(std_out,*) "compute_band_energy totch3",totch3
599 :
600 292 : if (occ_type == " lda") then
601 93 : if (abs(energies_dmft%eband_dft-energies_dmft%eband_dmft) > tol5) then
602 89 : write(message,'(5x,3a,15x,a,f12.6,a,15x,a,5x,f12.5)') "Warning: ", &
603 89 : & "Differences between band energy with Fermi-Dirac occupations",ch10, &
604 89 : & "and occupations from DFT Green's function is:",energies_dmft%eband_dft-energies_dmft%eband_dmft,ch10, &
605 178 : & "which is larger than",tol5
606 89 : call wrtout(std_out,message,'COLL')
607 : write(message,'(a)') &
608 89 : & " Action: increase the number of frequencies, or reduce the number of high energy DMFT bands"
609 89 : call wrtout(std_out,message,'COLL')
610 : else
611 4 : write(message,'(3a,10x,a,f12.6,a,10x,a,5x,f12.5)') " ", &
612 4 : & "Differences between band energy with Fermi-Dirac occupations",ch10, &
613 4 : & "and occupations from DFT Green's function is:",energies_dmft%eband_dft-energies_dmft%eband_dmft,ch10, &
614 8 : & "which is smaller than",tol5
615 4 : call wrtout(std_out,message,'COLL')
616 : end if ! tol
617 : end if ! occ_type=lda
618 :
619 292 : if (present(fcalc_dft)) then
620 0 : if (abs(energies_dmft%fband_dft-green%trace_log) > tol5) then
621 0 : write(message,'(5x,3a,15x,a,f12.6,a,15x,a,5x,f12.5)') "Warning: ", &
622 0 : & "Differences between free energy with Fermi-Dirac occupations",ch10, &
623 0 : & "and occupations from DFT Green's function is:",energies_dmft%fband_dft-green%trace_log,ch10, &
624 0 : & "which is larger than",tol5
625 0 : call wrtout(std_out,message,'COLL')
626 : write(message,'(a)') &
627 0 : & " Action: increase the number of frequencies, or reduce the number of high energy DMFT bands"
628 0 : call wrtout(std_out,message,'COLL')
629 : else
630 0 : write(message,'(3a,10x,a,f12.6,a,10x,a,5x,f12.5)') " ", &
631 0 : & "Differences between free energy with Fermi-Dirac occupations",ch10, &
632 0 : & "and occupations from DFT Green's function is:",energies_dmft%fband_dft-green%trace_log,ch10, &
633 0 : & "which is smaller than",tol5
634 0 : call wrtout(std_out,message,'COLL')
635 : end if ! tol
636 : end if ! occ_type=lda
637 :
638 292 : end subroutine compute_band_energy
639 : !!***
640 :
641 : !!****f* m_energy/compute_migdal_energy
642 : !! NAME
643 : !! compute_migdal_energy
644 : !!
645 : !! FUNCTION
646 : !! Computes Midgal energy = 1/2 * Tr(Sigma*G)
647 : !!
648 : !! INPUTS
649 : !! green <type(green_type)>= green function data
650 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
651 : !! self <type(self_type)>= self energy function data
652 : !! iatom = if present, only computes the contribution of this atom
653 : !!
654 : !! OUTPUT
655 : !! e_hu_migdal(natom)= Migdal energy for each atom.
656 : !! e_hu_mig_tot= Total Migdal energy.
657 : !!
658 : !! SOURCE
659 :
660 106 : subroutine compute_migdal_energy(e_hu_migdal,e_hu_migdal_tot,green,paw_dmft,self,iatom)
661 :
662 : !#ifdef FC_INTEL
663 : !DEC$ NOOPTIMIZE
664 : !#endif
665 :
666 : !Arguments ------------------------------------
667 : type(green_type), intent(in) :: green
668 : type(paw_dmft_type), intent(in) :: paw_dmft
669 : real(dp), intent(out) :: e_hu_migdal_tot
670 : real(dp), intent(inout) :: e_hu_migdal(paw_dmft%natom)
671 : type(self_type), target, intent(in) :: self
672 : integer, optional, intent(in) :: iatom
673 : ! integer :: prtopt
674 : !Local variables-------------------------------
675 : integer :: i,iatom_,ierr,ifreq,j,myproc,natom,nmoments,nspinor,nsppol,nwlo
676 : real(dp) :: beta,temp
677 : complex(dp) :: omega
678 106 : complex(dp), allocatable :: omega_fac(:),trace_moments(:,:),trace(:)
679 106 : type(matlu_type), allocatable :: self_nwlo_re(:)
680 106 : type(matlu_type), pointer :: matlu_tmp(:) => null()
681 : character(len=500) :: message
682 : ! *********************************************************************
683 :
684 : ! Only imaginary frequencies here
685 106 : if (green%w_type == "real" .or. self%w_type == "real") then
686 0 : message = 'compute_migdal_energy not implemented for real frequency'
687 0 : ABI_BUG(message)
688 : end if
689 :
690 : ! == Compute Correlation energy from Migdal formula
691 : ! -----------------------------------------------------------------------
692 106 : beta = one / paw_dmft%temp
693 106 : myproc = paw_dmft%myproc
694 106 : natom = paw_dmft%natom
695 106 : nmoments = 0
696 106 : nspinor = paw_dmft%nspinor
697 106 : nsppol = paw_dmft%nsppol
698 106 : nwlo = green%nw
699 106 : temp = paw_dmft%temp
700 :
701 106 : if (self%has_moments == 1) nmoments = self%nmoments
702 106 : iatom_ = 0
703 106 : if (present(iatom)) then
704 0 : iatom_ = iatom
705 0 : if (self%has_moments == 0) ABI_BUG("You should not be here")
706 : end if
707 :
708 106 : if (green%nw /= self%nw) then
709 0 : message = 'self and green do not contain the same number of frequencies'
710 0 : ABI_BUG(message)
711 : end if
712 : ! write(std_out,*) "beta",beta
713 :
714 424 : ABI_MALLOC(trace_moments,(natom,nmoments))
715 318 : ABI_MALLOC(trace,(natom))
716 :
717 476 : e_hu_migdal(:) = zero
718 476 : trace(:) = czero
719 :
720 106 : if (self%has_moments == 1) then
721 0 : trace_moments(:,:) = czero
722 0 : do i=1,nmoments
723 0 : do j=1,i
724 0 : call trace_prod_matlu(self%moments(j)%matlu(:),green%moments(i-j+1)%matlu(:),natom,trace(:),iatom=iatom_)
725 0 : trace_moments(:,i) = trace_moments(:,i) + trace(:)
726 : end do ! j
727 : end do ! i
728 : else
729 688 : ABI_MALLOC(self_nwlo_re,(natom))
730 582 : ABI_MALLOC(matlu_tmp,(natom))
731 106 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),matlu_tmp(:))
732 106 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),self_nwlo_re(:))
733 106 : call copy_matlu(self%oper(nwlo)%matlu(:),self_nwlo_re(:),natom,opt_re=1)
734 : end if ! moments
735 :
736 23195 : do ifreq=1,nwlo
737 :
738 23089 : if (self%distrib%procf(ifreq) /= myproc) cycle
739 :
740 8032 : omega = cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
741 :
742 8032 : if (self%has_moments == 1) then
743 0 : matlu_tmp => self%oper(ifreq)%matlu(:)
744 : else
745 8032 : call add_matlu(self%oper(ifreq)%matlu(:),self_nwlo_re(:),matlu_tmp(:),natom,-1)
746 : end if ! moments
747 :
748 8032 : call trace_prod_matlu(matlu_tmp(:),green%oper(ifreq)%matlu(:),natom,trace(:),iatom=iatom_)
749 :
750 38859 : e_hu_migdal(:) = e_hu_migdal(:) + dble(trace(:))*paw_dmft%wgt_wlo(ifreq)*temp*two
751 :
752 : end do ! ifreq
753 :
754 106 : call xmpi_sum(e_hu_migdal(:),paw_dmft%spacecomm,ierr)
755 :
756 318 : ABI_MALLOC(omega_fac,(nmoments))
757 :
758 106 : do i=1,nmoments
759 0 : omega_fac(i) = czero
760 0 : do ifreq=nwlo,1,-1 ! NEVER change the summation order and DON'T use the intrinsic SUM
761 0 : omega_fac(i) = omega_fac(i) + cone / (paw_dmft%omega_lo(ifreq))**i
762 : end do
763 0 : omega_fac(i) = - two * temp * omega_fac(i) / (j_dpc)**i
764 0 : if (i == 1) omega_fac(i) = omega_fac(i) + half
765 0 : if (i == 2) omega_fac(i) = omega_fac(i) - cone/(four*temp)
766 0 : if (i == 4) omega_fac(i) = omega_fac(i) + cone/(dble(48)*(temp**3))
767 106 : e_hu_migdal(:) = e_hu_migdal(:) + dble(trace_moments(:,i)*omega_fac(i))
768 : end do
769 :
770 106 : if (self%has_moments /= 1) then
771 106 : call destroy_matlu(matlu_tmp(:),natom)
772 106 : ABI_FREE(matlu_tmp)
773 : end if
774 106 : matlu_tmp => null()
775 :
776 106 : ABI_FREE(omega_fac)
777 :
778 106 : if (self%has_moments == 0) then
779 106 : call trace_prod_matlu(self_nwlo_re(:),green%occup%matlu(:),natom,trace(:))
780 106 : call destroy_matlu(self_nwlo_re(:),natom)
781 476 : ABI_FREE(self_nwlo_re)
782 476 : e_hu_migdal(:) = e_hu_migdal(:) + dble(trace(:))
783 : end if
784 :
785 106 : ABI_FREE(trace_moments)
786 :
787 476 : e_hu_migdal(:) = half * e_hu_migdal(:) ! E_mig = 1/2 * Tr(Sig*G)
788 476 : e_hu_migdal_tot = sum(e_hu_migdal(:))
789 :
790 106 : ABI_FREE(trace)
791 :
792 : !xmig_1=zero
793 : !xmig_2=zero
794 : !xmig_3=zero
795 :
796 : !e_hu_migdal_tot = zero
797 : !do iatom=1,natom
798 : ! shift=czero
799 : ! if(paw_dmft%dmft_solv==4) shift=self%qmc_shift(iatom)+self%qmc_xmu(iatom)
800 : ! write(std_out,*) "shiftttt",shift
801 : ! lpawu=paw_dmft%lpawu(iatom)
802 : ! if(lpawu/=-1) then
803 : ! xmig_1=czero
804 : ! xmig_2=czero
805 : ! xmig_3=czero
806 : ! ndim=2*lpawu+1
807 : ! do isppol=1,nsppol
808 : ! do ispinor = 1 , nspinor
809 : ! do ispinor1 = 1, nspinor
810 : ! do im=1,ndim
811 : ! do im1=1,ndim
812 : ! do ifreq=1,nwlo
813 : ! write(std_out,*) ifreq,xmig_1,imag(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)),&
814 : !& green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )
815 : ! xmig_1=xmig_1 + j_dpc/beta* &
816 : !& aimag(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))* &
817 : !& green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )* &
818 : !& paw_dmft%wgt_wlo(ifreq)
819 : ! if(ispinor==ispinor1.and.im==im1) then
820 : ! se=(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)- &
821 : !& self%oper (nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))
822 : ! else
823 : ! se=self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)
824 : ! endif
825 : ! xmig_2=xmig_2 + one/beta*real(se)* &
826 : !& green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )* &
827 : !& paw_dmft%wgt_wlo(ifreq)
828 : ! if(ispinor==ispinor1.nd.im==im1.and.ifreq==1) then
829 : ! if(ifreq==1) then
830 : ! xmig_3=xmig_3 + &
831 : !& real(self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)+shift)* &
832 : !& green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two
833 : ! write(std_out,*) "xmig_3",xmig_3
834 : ! write(std_out,*) "self",self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)
835 : ! write(std_out,*) "shift",shift
836 : ! write(std_out,*) "occup", green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two
837 : ! endif
838 : ! enddo
839 : ! if(ispinor==ispinor1.and.im==im1) then
840 : ! xmig_3=xmig_3 + &
841 : !& real(self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))* &
842 : !!& green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two
843 : ! endif
844 : ! enddo
845 : ! enddo
846 : ! enddo
847 : ! enddo
848 : ! enddo
849 : ! if(nsppol==1.and.nspinor==1) then
850 : ! e_hu_migdal(iatom)=two*real(xmig_1+xmig_2+xmig_3)
851 : ! else
852 : ! e_hu_migdal(iatom)=real(xmig_1+xmig_2+xmig_3)
853 : ! endif
854 : ! e_hu_migdal_tot = e_hu_migdal_tot + e_hu_migdal(iatom)
855 : ! if(abs(pawprtvol)>=3) then
856 : ! write(message,'(2a,3(a,5x,a,2f12.6))')ch10,&
857 : !& " Interaction energy: Decomposition of Migdal energy",ch10,&
858 : !& "xmig_1=",xmig_1,ch10,&
859 : !& "xmig_2=",xmig_2,ch10,&
860 : !& "xmig_3=",xmig_3
861 : ! call wrtout(std_out,message,'COLL')
862 : ! endif
863 : ! endif ! lpawu
864 : !enddo
865 :
866 106 : end subroutine compute_migdal_energy
867 : !!***
868 :
869 : !!****f* m_energy/compute_dftu_energy
870 : !! NAME
871 : !! compute_dftu_energy
872 : !!
873 : !! FUNCTION
874 : !! Initialize noccmmp from green%occup and compute DFT+U energy with it
875 : !!
876 : !! INPUTS
877 : !! energies_dmft = datastructure for dmft energy
878 : !! green <type(green_type)>= green function data
879 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
880 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
881 : !! renorm = if present change U->1 and J-> renorm just for pawuenergy
882 : !! renorm = J/U for the real values (does not depend on "lambda" entropy)
883 : !!
884 : !! OUTPUT
885 : !!
886 : !! SOURCE
887 :
888 208 : subroutine compute_dftu_energy(energies_dmft,green,paw_dmft,pawtab,renorm)
889 :
890 : !Arguments ------------------------------------
891 : !type
892 : type(energy_type), intent(inout) :: energies_dmft
893 : type(green_type), intent(in) :: green
894 : type(paw_dmft_type), intent(in) :: paw_dmft
895 : type(pawtab_type), intent(in) :: pawtab(paw_dmft%ntypat)
896 : real(dp), optional, intent(in) :: renorm(:)
897 : ! integer :: prtopt
898 : !Local variables-------------------------------
899 : integer :: iatom,idijeff,im,im1,ims,ims1,ispinor,ispinor1,isppol,itypat
900 : integer :: lpawu,lpawu1,ndim,ndim1,nocc,nsploop,prt_pawuenergy
901 : real(dp) :: e_dc,e_dc_for_s,e_dcdc,e_dcdc_for_s,e_ee,edftumdc,edftumdc_for_s
902 : real(dp) :: edftumdcdc,edftumdcdc_for_s,e_ee_for_s,jpawu,upawu,xe1,xe2
903 : logical :: dmft_optim,t2g,x2my2d
904 : character(len=500) :: message
905 : integer, parameter :: spinor_idxs(2,4) = RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
906 : integer, parameter :: mt2g(3) = (/1,2,4/)
907 208 : real(dp), allocatable :: noccmmp(:,:,:,:),nocctot(:)
908 : ! *********************************************************************
909 :
910 : ! - allocations
911 : ! -----------------------------------------------------------------------
912 :
913 208 : e_dc = zero
914 208 : e_dc_for_s = zero
915 208 : e_dcdc = zero
916 208 : e_ee = zero
917 208 : edftumdc = zero
918 208 : edftumdcdc = zero
919 208 : nsploop = max(paw_dmft%nsppol,paw_dmft%nspinor**2)
920 208 : nocc = nsploop
921 208 : t2g = (paw_dmft%dmft_t2g == 1)
922 208 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
923 :
924 208 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
925 :
926 208 : isppol = 1
927 208 : ispinor = 1
928 208 : ispinor1 = 1
929 :
930 624 : ABI_MALLOC(nocctot,(nocc))
931 :
932 : ! - Loop and call to pawuenergy
933 : ! -----------------------------------------------------------------------
934 932 : do iatom=1,paw_dmft%natom
935 724 : lpawu = paw_dmft%lpawu(iatom)
936 724 : if (lpawu == -1) cycle
937 244 : itypat = paw_dmft%typat(iatom)
938 244 : lpawu1 = lpawu
939 244 : if (t2g .or. x2my2d) lpawu1 = 2
940 244 : ndim = 2*lpawu + 1
941 244 : ndim1 = 2*lpawu1 + 1
942 :
943 1220 : ABI_MALLOC(noccmmp,(2,ndim1,ndim1,nocc))
944 35041 : noccmmp(:,:,:,:) = zero
945 :
946 : ! - Setup nocctot and noccmmp
947 : ! -----------------------------------------------------------------------
948 627 : nocctot(:) = zero ! contains nmmp in the n m representation
949 : ! Begin loop over spin/spinors to initialize noccmmp
950 627 : do idijeff=1,nsploop
951 :
952 383 : if (nsploop <= 2) then
953 : isppol = idijeff
954 12 : else if (nsploop == 4) then
955 12 : ispinor = spinor_idxs(1,idijeff)
956 12 : ispinor1 = spinor_idxs(2,idijeff)
957 : else
958 0 : write(message,'(2a)') " BUG in m_energy: nsploop should be equal to 1, 2 or 4"
959 0 : call wrtout(std_out,message,'COLL')
960 : end if ! nsploop
961 : ! Initialize noccmmp
962 2160 : do im1=1,ndim
963 1777 : ims1 = im1
964 : ! Correct bug in computation of DFT+U energy in the t2g/x2my2d case with TRIQS
965 1777 : if (x2my2d .and. dmft_optim) ims1 = 5
966 1777 : if (t2g .and. dmft_optim) ims1 = mt2g(im1)
967 11039 : do im=1,ndim
968 8879 : ims = im
969 8879 : if (x2my2d .and. dmft_optim) ims = 5
970 8879 : if (t2g .and. dmft_optim) ims = mt2g(im)
971 : ! Here, we take the transpose in order to match pawuenergy's conventions
972 : noccmmp(1,ims,ims1,idijeff) = &
973 8879 : & dble(green%occup%matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol))
974 : noccmmp(2,ims,ims1,idijeff) = &
975 10656 : & aimag(green%occup%matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol))
976 : ! noccmmp(1,im,im1,idijeff)=real(green%occup%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
977 : ! noccmmp(2,im,im1,idijeff)=imag(green%occup%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
978 : end do ! im
979 : end do ! im1
980 : ! Compute nocctot
981 627 : if (green%has_charge_matlu_solver /= 2) then
982 179 : if (nsploop == 4 .and. idijeff <= 2) then
983 16 : do im1=1,ndim
984 : nocctot(1) = nocctot(1) + &
985 16 : & dble(green%occup%matlu(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol))
986 : end do ! im1
987 177 : else if (nsploop <= 2) then
988 960 : do im1=1,ndim
989 : nocctot(idijeff) = nocctot(idijeff) + &
990 960 : & dble(green%occup%matlu(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol))
991 : end do ! im1
992 : end if ! nsploop
993 : else
994 204 : if (nsploop == 4 .and. idijeff == 1) then
995 2 : nocctot(1) = green%charge_matlu_solver(2,iatom) ! total nb of elec for nspinor=2 is (iatom,2) !!
996 202 : else if (nsploop <= 2) then
997 196 : nocctot(idijeff) = green%charge_matlu_solver(idijeff,iatom) ! first spin
998 : end if ! nsploop
999 : end if ! charge_matlu_solver
1000 : end do ! idijeff
1001 :
1002 244 : xe1 = e_dc
1003 244 : xe2 = e_ee
1004 : ! write(std_out,*)" nocctot(1)",nocctot(1),green%charge_matlu_solver(iatom,1)
1005 244 : edftumdc = zero
1006 244 : edftumdcdc = zero
1007 244 : if (present(renorm)) then
1008 9 : upawu = one
1009 9 : jpawu = renorm(iatom)
1010 9 : prt_pawuenergy = 0
1011 : else
1012 235 : upawu = pawtab(itypat)%upawu
1013 235 : jpawu = pawtab(itypat)%jpawu
1014 235 : prt_pawuenergy = 3
1015 : end if ! present(renorm)
1016 :
1017 : call pawuenergy(iatom,edftumdc,edftumdcdc,noccmmp(:,:,:,:),nocctot(:),prt_pawuenergy,pawtab(itypat),dmft_dc=paw_dmft%dmft_dc,&
1018 244 : & e_ee=e_ee,e_dc=e_dc,e_dcdc=e_dcdc,u_dmft=upawu,j_dmft=jpawu,paw_dmft=paw_dmft)
1019 :
1020 244 : if (paw_dmft%ientropy == 1) then
1021 : call pawuenergy(iatom,edftumdc_for_s,edftumdcdc_for_s,noccmmp(:,:,:,:),nocctot(:),prt_pawuenergy, &
1022 : & pawtab(itypat),dmft_dc=paw_dmft%dmft_dc,e_ee=e_ee_for_s,e_dc=e_dc_for_s,e_dcdc=e_dcdc_for_s,&
1023 0 : & u_dmft=paw_dmft%u_for_s/Ha_eV,j_dmft=paw_dmft%j_for_s/Ha_eV)
1024 : end if
1025 :
1026 244 : energies_dmft%e_dc(iatom) = e_dc - xe1
1027 244 : energies_dmft%e_hu_dftu(iatom) = e_ee - xe2
1028 :
1029 932 : ABI_FREE(noccmmp)
1030 : end do ! iatom
1031 :
1032 208 : ABI_FREE(nocctot)
1033 :
1034 : ! - gather results
1035 : ! -----------------------------------------------------------------------
1036 208 : energies_dmft%e_dc_tot = e_dc ! this is the only quantity used afterwards.
1037 208 : energies_dmft%e_hu_dftu_tot = e_ee
1038 208 : energies_dmft%e_dcdc = e_dcdc
1039 208 : if (paw_dmft%ientropy == 1) then
1040 0 : write(message,'(a,3(f14.10,3x))') "For entropy calculation E_dc_tot, u_for_s, j_for,s", &
1041 0 : & e_dc_for_s,paw_dmft%u_for_s,paw_dmft%j_for_s
1042 0 : call wrtout(std_out,message,'COLL')
1043 0 : write(message,'(a,3(f14.10,3x))') "Reference calculation E_dc_tot, upawu , jpawu ",&
1044 0 : & e_dc,upawu*Ha_eV,jpawu*Ha_eV
1045 0 : call wrtout(std_out,message,'COLL')
1046 : end if ! ientropy=1
1047 :
1048 208 : end subroutine compute_dftu_energy
1049 : !!***
1050 :
1051 : !!****f* m_energy/compute_noninterentropy
1052 : !! NAME
1053 : !! compute_noninterentropy
1054 : !!
1055 : !! FUNCTION
1056 : !!
1057 : !! INPUTS
1058 : !! cryst_struc <type(crystal_t)>=crystal structure data
1059 : !! green <type(green_type)>= green function data only for Tr(G(self-hdc))
1060 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1061 : !!
1062 : !! OUTPUT
1063 : !!
1064 : !! SIDE EFFECTS
1065 : !!
1066 : !! SOURCE
1067 :
1068 0 : subroutine compute_noninterentropy(cryst_struc,green,paw_dmft)
1069 :
1070 : use m_crystal, only : crystal_t
1071 :
1072 : !Arguments ------------------------------------
1073 : !type
1074 : type(crystal_t),intent(in) :: cryst_struc
1075 : type(green_type),intent(in) :: green
1076 : type(paw_dmft_type), intent(inout) :: paw_dmft
1077 :
1078 : !Local variables-------------------------------
1079 : integer :: ib,ikpt,isppol,natom,nspinor,nsppol
1080 : real(dp) :: beta,eig,fermi,s_1,s_2,occ1,occ2,f_1,e_1,f_1a,s_1a,e_2
1081 : character(len=800) :: message
1082 : ! *********************************************************************
1083 0 : write(message,'(2a,i6)') ch10," == Compute T*Entropy for fermi level and DFT-KS eigenvalues "
1084 0 : call wrtout(std_out,message,'COLL')
1085 :
1086 0 : natom=cryst_struc%natom
1087 0 : nsppol=paw_dmft%nsppol
1088 0 : nspinor=paw_dmft%nspinor
1089 0 : beta=one/paw_dmft%temp
1090 0 : s_1=zero
1091 0 : s_1a=zero
1092 0 : f_1=zero
1093 0 : f_1a=zero
1094 0 : e_1=zero
1095 0 : e_2=zero
1096 0 : s_2=zero
1097 0 : do isppol=1,paw_dmft%nsppol
1098 0 : do ikpt=1,paw_dmft%nkpt
1099 0 : do ib=1,paw_dmft%mbandc
1100 0 : eig=paw_dmft%eigen_dft(ib,ikpt,isppol)
1101 : fermi=paw_dmft%fermie_dft
1102 0 : fermi=paw_dmft%fermie
1103 0 : occ1=occup_fd(eig,fermi,paw_dmft%temp)
1104 0 : occ2=green%occup%ks(ib,ib,ikpt,isppol)
1105 : ! write(std_out,*) occ1,occ2
1106 :
1107 : ! entropy from Fermi Dirac
1108 0 : if((occ1.ge.tol9).and.((one-occ1).ge.tol9)) then
1109 0 : s_1=s_1+(occ1*log(occ1)+(one-occ1)*log(one-occ1))*paw_dmft%wtk(ikpt)
1110 : ! write(std_out,*) occ1,one-occ1,"p1",(occ1*log(occ1)+(one-occ1)*log(one-occ1))*paw_dmft%wtk(ikpt)
1111 : endif
1112 :
1113 : ! Free energy from Fermi Dirac
1114 0 : if((eig-fermi).ge.zero) then ! occ1 -> 0 ; 1-occ1 -> 1
1115 0 : f_1=f_1-paw_dmft%wtk(ikpt)/beta*log(one+exp(-beta*(eig-fermi)))
1116 0 : f_1a=f_1a+paw_dmft%wtk(ikpt)/beta*log(one-occ1)
1117 0 : s_1a=s_1a+((one-occ1)*log(one-occ1)+occ1*(-beta*(eig-fermi)+log(one-occ1)))*paw_dmft%wtk(ikpt)
1118 : else ! occ1 -> 1 , 1-occ1 -> 0
1119 0 : f_1=f_1-paw_dmft%wtk(ikpt)/beta*(log(one+exp(beta*(eig-fermi)))-beta*(eig-fermi))
1120 0 : f_1a=f_1a-paw_dmft%wtk(ikpt)/beta*(-log(occ1)-beta*(eig-fermi))
1121 0 : s_1a=s_1a+(occ1*log(occ1)+(one-occ1)*(beta*(eig-fermi)+log(occ1)))*paw_dmft%wtk(ikpt)
1122 : endif
1123 :
1124 : ! Internal energy from Fermi Dirac
1125 0 : e_1=e_1+(eig-fermi)*paw_dmft%wtk(ikpt)*occ1
1126 0 : e_2=e_2+(eig-fermi)*paw_dmft%wtk(ikpt)*occ2
1127 :
1128 : ! entropy from green function occupations.
1129 0 : if((occ2.ge.tol9).and.((one-occ2).ge.tol9)) then
1130 : ! write(std_out,*) occ2,one-occ2,"p2",(occ2*log(occ2)+(one-occ2)*log(one-occ2))*paw_dmft%wtk(ikpt)
1131 0 : s_2=s_2+(occ2*log(occ2)+(one-occ2)*log(one-occ2))*paw_dmft%wtk(ikpt)
1132 : endif
1133 : enddo
1134 : enddo
1135 : enddo
1136 0 : s_1=-s_1*paw_dmft%temp
1137 0 : s_1a=-s_1a*paw_dmft%temp
1138 0 : s_2=-s_2*paw_dmft%temp
1139 :
1140 : write(message,'(8(2a,e20.9))') &
1141 0 : & ch10," T*Entropy from Fermi Dirac occupations ", s_1,&
1142 0 : & ch10," T*Entropy from Fermi Dirac occupations 2 ", s_1a,&
1143 0 : & ch10," T*Entropy from Green function occupations ", s_2,&
1144 0 : & ch10," Free energy F ", f_1,&
1145 0 : & ch10," Free energy Fa ", f_1a,&
1146 0 : & ch10," internal energy U ", e_1,&
1147 0 : & ch10," internal energy U from Gr Func Occ ", e_2,&
1148 0 : & ch10," U-F ", e_1-f_1
1149 0 : call wrtout(std_out,message,'COLL')
1150 :
1151 :
1152 0 : end subroutine compute_noninterentropy
1153 : !!***
1154 :
1155 : !!****f* m_energy/compute_free_energy
1156 : !! NAME
1157 : !! compute_free_energy
1158 : !!
1159 : !! FUNCTION
1160 : !! Computes the different DFT+DMFT contributions to the free energy.
1161 : !!
1162 : !! INPUTS
1163 : !! energies_dmft = datastructure for dmft energy
1164 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1165 : !! green <type(green_type)>= green function data
1166 : !! self <type(self_type)>= self energy function data
1167 : !! weiss <type(green_type)>= weiss function data
1168 : !! part = "band" : computes Tr(log(G_DFT)) in KS space
1169 : !! = "main" : computes Tr(Sig*G) and Tr(log(G))
1170 : !! = "impu" : computes the rest
1171 : !!
1172 : !! OUTPUT
1173 : !!
1174 : !! SIDE EFFECTS
1175 : !!
1176 : !! SOURCE
1177 :
1178 0 : subroutine compute_free_energy(energies_dmft,paw_dmft,green,part,self)
1179 :
1180 : !Arguments ------------------------------------
1181 : type(energy_type), intent(inout) :: energies_dmft
1182 : type(paw_dmft_type), intent(in) :: paw_dmft
1183 : type(green_type), intent(in) :: green
1184 : type(self_type), optional, intent(in) :: self
1185 : character(len=4), intent(in) :: part
1186 : !Local variables-------------------------------
1187 : integer :: integral
1188 0 : real(dp), allocatable :: e_hu_tmp(:)
1189 : ! *********************************************************************
1190 :
1191 0 : ABI_MALLOC(e_hu_tmp,(paw_dmft%natom))
1192 :
1193 0 : integral = paw_dmft%dmft_triqs_compute_integral
1194 :
1195 : ! Compare Tr(log(G_DFT)) with the analytical formula (not used, simply to check that we have enough frequencies)
1196 0 : if (part == "band") then
1197 0 : call compute_band_energy(energies_dmft,green,paw_dmft,"nlda",fcalc_dft=1)
1198 : end if
1199 :
1200 0 : if (part == "impu") then
1201 :
1202 : ! Ekin_imp
1203 0 : energies_dmft%ekin_imp = green%ekin_imp
1204 :
1205 0 : if (integral == 1) then
1206 0 : energies_dmft%fband_weiss = green%fband_weiss
1207 0 : energies_dmft%fimp = energies_dmft%fband_weiss + green%integral
1208 : else
1209 0 : energies_dmft%fimp = energies_dmft%ekin_imp + energies_dmft%e_hu_tot
1210 : end if ! integral
1211 :
1212 : ! Integral of <dH/dlambda>
1213 0 : if (integral > 0) energies_dmft%integral = green%integral
1214 :
1215 : ! Tr(log(G_imp))
1216 0 : call compute_trace_log_loc(green,paw_dmft,energies_dmft%fband_imp)
1217 :
1218 : ! Tr(Sigma_imp*G_imp)
1219 0 : call compute_migdal_energy(e_hu_tmp(:),energies_dmft%emig_imp,green,paw_dmft,self)
1220 0 : energies_dmft%emig_imp = two * energies_dmft%emig_imp
1221 :
1222 : ! simp = entropy of the impurity
1223 0 : energies_dmft%simp = (energies_dmft%ekin_imp+energies_dmft%e_hu_tot-energies_dmft%fimp) / paw_dmft%temp
1224 :
1225 : end if ! part="impu"
1226 :
1227 0 : if (part == "main") then
1228 :
1229 : ! Tr(Sigma*G)
1230 0 : call compute_migdal_energy(e_hu_tmp(:),energies_dmft%emig_loc,green,paw_dmft,self)
1231 0 : energies_dmft%emig_loc = two * energies_dmft%emig_loc
1232 :
1233 : ! Tr(log(G)) + mu*N
1234 0 : energies_dmft%fband_dmft = green%trace_log
1235 :
1236 : ! fdmft = F_{dft+dmft} - E_dft
1237 : energies_dmft%fdmft = energies_dmft%fband_dmft - energies_dmft%eband_dmft - energies_dmft%emig_loc &
1238 0 : & - energies_dmft%fband_imp + energies_dmft%emig_imp + energies_dmft%fimp - energies_dmft%e_dcdc
1239 :
1240 : ! sdmft = total entropy
1241 0 : energies_dmft%sdmft = (energies_dmft%edmft-energies_dmft%fdmft) / paw_dmft%temp
1242 :
1243 0 : call print_free_energy(energies_dmft,paw_dmft)
1244 :
1245 : end if ! part="main"
1246 :
1247 0 : ABI_FREE(e_hu_tmp)
1248 :
1249 0 : end subroutine compute_free_energy
1250 : !!***
1251 :
1252 : !!****f* m_energy/compute_trace_log_loc
1253 : !! NAME
1254 : !! compute_trace_log_loc
1255 : !!
1256 : !! FUNCTION
1257 : !! Computes Tr(log(G)) in local space.
1258 : !!
1259 : !! INPUTS
1260 : !! green <type(green_type)>= green function data
1261 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1262 : !! opt_inv = 0 (default) when green = G
1263 : !! = 1 when green = G^-1 (used for Weiss field)
1264 : !! CAREFUL: in this case, the chemical potential is shifted
1265 : !! by dmft_triqs_shift_mu
1266 : !!
1267 : !! OUTPUT
1268 : !! trace = Tr(log(G_loc))
1269 : !!
1270 : !! SIDE EFFECTS
1271 : !!
1272 : !! SOURCE
1273 :
1274 0 : subroutine compute_trace_log_loc(green,paw_dmft,trace,opt_inv)
1275 :
1276 : !Arguments ------------------------------------
1277 : type(green_type), target, intent(in) :: green
1278 : type(paw_dmft_type), intent(in) :: paw_dmft
1279 : real(dp), intent(out) :: trace
1280 : integer, optional, intent(in) :: opt_inv
1281 : !Local variables-------------------------------
1282 : integer :: i,iatom,ierr,ifreq,im,info,isppol,lpawu,lwork,natom,ndim
1283 : integer :: nmoments,nspinor,nsppol,nwlo,optinv
1284 : real(dp) :: correction,fac,freq2,temp
1285 : complex(dp) :: trace_tmp
1286 0 : real(dp), allocatable :: eig(:),rwork(:)
1287 0 : complex(dp), allocatable :: mat_temp(:,:),omega_fac(:),work(:)
1288 0 : complex(dp), target, allocatable :: mat_temp2(:,:)
1289 : complex(dp), pointer :: mat_pt(:,:) => null()
1290 : ! *********************************************************************
1291 :
1292 0 : optinv = 0
1293 0 : if (present(opt_inv)) optinv = opt_inv
1294 :
1295 0 : natom = paw_dmft%natom
1296 0 : nmoments = green%nmoments - 1
1297 0 : nspinor = paw_dmft%nspinor
1298 0 : nsppol = paw_dmft%nsppol
1299 0 : nwlo = green%nw
1300 0 : temp = paw_dmft%temp
1301 0 : trace = zero
1302 0 : ndim = nspinor * (2*paw_dmft%maxlpawu+1)
1303 :
1304 0 : ABI_MALLOC(eig,(ndim))
1305 0 : ABI_MALLOC(rwork,(3*ndim-2))
1306 0 : ABI_MALLOC(work,(2*ndim-1))
1307 0 : ABI_MALLOC(mat_temp,(ndim,ndim))
1308 0 : call zheev('n','u',ndim,mat_temp(:,:),ndim,eig(:),work(:),-1,rwork(:),info)
1309 0 : lwork = int(work(1))
1310 0 : ABI_FREE(work)
1311 0 : ABI_MALLOC(work,(lwork))
1312 0 : ABI_FREE(mat_temp)
1313 :
1314 0 : do ifreq=1,nwlo
1315 0 : if (green%distrib%procf(ifreq) /= paw_dmft%myproc) cycle
1316 0 : fac = merge(temp*two,temp,nsppol==1.and.nspinor==1)
1317 0 : trace_tmp = czero
1318 0 : freq2 = paw_dmft%omega_lo(ifreq)**2
1319 0 : do iatom=1,natom
1320 0 : lpawu = paw_dmft%lpawu(iatom)
1321 0 : if (lpawu == -1) cycle
1322 0 : ndim = nspinor * (2*lpawu+1)
1323 0 : ABI_MALLOC(mat_temp,(ndim,ndim))
1324 0 : do isppol=1,nsppol
1325 :
1326 0 : if (optinv == 0) then
1327 0 : mat_pt => green%oper(ifreq)%matlu(iatom)%mat(:,:,isppol)
1328 0 : else if (optinv == 1) then
1329 0 : ABI_MALLOC(mat_temp2,(ndim,ndim))
1330 0 : mat_temp2(:,:) = green%oper(ifreq)%matlu(iatom)%mat(:,:,isppol)
1331 0 : do im=1,ndim
1332 0 : mat_temp2(im,im) = mat_temp2(im,im) + paw_dmft%dmft_triqs_shift_mu
1333 : end do ! im
1334 0 : mat_pt => mat_temp2(:,:)
1335 : end if
1336 :
1337 : call abi_xgemm("n","c",ndim,ndim,ndim,cone,mat_pt(:,:),ndim,mat_pt(:,:),ndim, &
1338 0 : & czero,mat_temp(:,:),ndim)
1339 0 : call zheev('n','u',ndim,mat_temp(:,:),ndim,eig(:),work(:),lwork,rwork(1:3*ndim-2),info)
1340 :
1341 0 : if (optinv == 1) then
1342 0 : trace_tmp = trace_tmp - sum(log(eig(1:ndim)/freq2))
1343 : else
1344 0 : trace_tmp = trace_tmp + sum(log(eig(1:ndim)*freq2))
1345 : end if
1346 :
1347 0 : ABI_SFREE(mat_temp2)
1348 :
1349 : end do ! isppol
1350 0 : if (ifreq == nwlo) then
1351 0 : correction = fac * nsppol * ndim * log(two)
1352 0 : trace = trace - correction
1353 : end if
1354 0 : ABI_FREE(mat_temp)
1355 : end do ! iatom
1356 0 : trace = trace + dble(trace_tmp)*fac
1357 : end do ! ifreq
1358 :
1359 0 : mat_pt => null()
1360 0 : ABI_FREE(rwork)
1361 0 : ABI_FREE(work)
1362 0 : ABI_FREE(eig)
1363 :
1364 0 : call xmpi_sum(trace,paw_dmft%spacecomm,ierr)
1365 :
1366 0 : ABI_MALLOC(omega_fac,(nmoments))
1367 :
1368 0 : do i=1,nmoments
1369 0 : omega_fac(i) = czero
1370 0 : do ifreq=nwlo,1,-1 ! NEVER change the summation order and DON'T use the intrinsic SUM
1371 0 : omega_fac(i) = omega_fac(i) + cone / (paw_dmft%omega_lo(ifreq))**i
1372 : end do
1373 0 : omega_fac(i) = - two * temp * omega_fac(i) / (j_dpc)**i
1374 0 : if (i == 1) omega_fac(i) = omega_fac(i) + half
1375 0 : if (i == 2) omega_fac(i) = omega_fac(i) - cone/(four*temp)
1376 0 : if (i == 4) omega_fac(i) = omega_fac(i) + cone/(dble(48)*(temp**3))
1377 : end do ! i
1378 :
1379 : ! Do not use dot_product
1380 0 : trace = trace + dble(sum(green%trace_moments_log_loc(1:nmoments)*omega_fac(1:nmoments)))
1381 :
1382 0 : ABI_FREE(omega_fac)
1383 :
1384 0 : end subroutine compute_trace_log_loc
1385 : !!***
1386 :
1387 : !!****f* m_energy/print_free_energy
1388 : !! NAME
1389 : !! print_free_energy
1390 : !!
1391 : !! FUNCTION
1392 : !! Prints the different DFT+DMFT contributions to the free energy.
1393 : !!
1394 : !! INPUTS
1395 : !! energies_dmft = datastructure for dmft energy
1396 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1397 : !!
1398 : !! OUTPUT
1399 : !!
1400 : !! SIDE EFFECTS
1401 : !!
1402 : !! SOURCE
1403 :
1404 0 : subroutine print_free_energy(energies_dmft,paw_dmft)
1405 :
1406 : !Arguments ------------------------------------
1407 : type(energy_type), intent(in) :: energies_dmft
1408 : type(paw_dmft_type), intent(in) :: paw_dmft
1409 : !Local variables-------------------------------
1410 : integer :: integral
1411 : real(dp) :: temp
1412 : character(len=10000) :: message,message2
1413 : ! *********************************************************************
1414 :
1415 0 : integral = paw_dmft%dmft_triqs_compute_integral
1416 0 : temp = paw_dmft%temp
1417 :
1418 0 : write(message,'(a,5x,2a,5x,a,10(a,5x,a,2x,f18.11),a)') ch10, &
1419 0 : & "-----------------------------------------------",ch10, &
1420 0 : & "--- Free Energy in DMFT (in Ha) ",ch10, &
1421 0 : & "--- E_hu (1) (Ha.) = ",energies_dmft%e_hu_tot,ch10, &
1422 0 : & "--- E_dc (2) (Ha.) = ",energies_dmft%e_dc_tot,ch10, &
1423 0 : & "--- E_dmft (1)-(2) (Ha.) = ",energies_dmft%edmft,ch10, &
1424 0 : & "--- Tr(log(G))+mu*N (3) (Ha.) = ",energies_dmft%fband_dmft,ch10, &
1425 0 : & "--- Tr(Sigma*G) (4) (Ha.) = ",energies_dmft%emig_loc,ch10, &
1426 0 : & "--- Tr(V_dc*rho) (5) (Ha.) = ",-energies_dmft%e_dcdc+energies_dmft%e_dc_tot,ch10, &
1427 0 : & "--- E_band_dmft (6) (Ha.) = ",energies_dmft%eband_dmft,ch10, &
1428 0 : & "--- Tr(log(Gimp)) (7) (Ha.) = ",energies_dmft%fband_imp,ch10, &
1429 0 : & "--- Tr(Sigma_imp*G_imp) (8) (Ha.) = ",energies_dmft%emig_imp,ch10, &
1430 0 : & "--- E_kinetic_imp (9) (Ha.) = ",energies_dmft%ekin_imp,ch10
1431 :
1432 0 : if (integral > 0) then
1433 : write(message2,'(8(5x,a,2x,f18.11,a),5x,a)') &
1434 0 : & "--- Tr(log(G0)) (10) (Ha.) = ",energies_dmft%fband_weiss,ch10, &
1435 0 : & "--- Integral (11) (Ha.) = ",energies_dmft%integral,ch10, &
1436 0 : & "--- F_imp (10)+(11) (Ha.) = ",energies_dmft%fimp,ch10, &
1437 0 : & "--- (-kT)*S_imp (10)+(11)-(9)-(1) (Ha.) = ",-temp*energies_dmft%simp,ch10, &
1438 0 : & "--- S_imp (12) = ",energies_dmft%simp,ch10, &
1439 0 : & "--- F_dmft (13) (Ha.) = ",energies_dmft%fdmft,ch10, &
1440 0 : & "--- (-kT)*S_dmft (13)-(1)+(2) (Ha.) = ",-temp*energies_dmft%sdmft,ch10, &
1441 0 : & "--- S_dmft = ",energies_dmft%sdmft,ch10, &
1442 0 : & "-----------------------------------------------"
1443 : else
1444 : write(message2,'(3(5x,a,2x,f18.11,a),5x,a)') &
1445 0 : & "--- F_dmft+T*S_imp (10) (Ha.) = ",energies_dmft%fdmft,ch10, &
1446 0 : & "--- (-kT)*(S_dmft-S_imp) (11) (Ha.) = ",-temp*energies_dmft%sdmft,ch10, &
1447 0 : & "--- S_dmft-S_imp (12) = ",energies_dmft%sdmft,ch10, &
1448 0 : & "-----------------------------------------------"
1449 : end if ! integral
1450 :
1451 0 : call wrtout(std_out,trim(adjustl(message))//trim(message2),'COLL')
1452 :
1453 0 : end subroutine print_free_energy
1454 : !!***
1455 :
1456 0 : END MODULE m_energy
1457 : !!***
|
