Branch data 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), 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 : 106 : 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 [ - + + - : 318 : ABI_MALLOC(energies_dmft%e_dc,(natom))
+ - + - -
+ - + -
+ ]
158 [ - + + - : 212 : ABI_MALLOC(energies_dmft%e_hu_dftu,(natom))
- + - + -
+ ]
159 [ - + + - : 212 : ABI_MALLOC(energies_dmft%e_hu_mig,(natom))
- + - + -
+ ]
160 [ - + + - : 212 : ABI_MALLOC(energies_dmft%e_hu_qmc,(natom))
- + - + -
+ ]
161 : 106 : energies_dmft%e_hu => energies_dmft%e_hu_mig(:)
162 [ + + ]: 476 : energies_dmft%e_dc(:) = zero
163 [ + + ]: 476 : energies_dmft%e_hu_dftu(:) = zero
164 [ + + ]: 476 : energies_dmft%e_hu_mig(:) = zero
165 [ + + ]: 476 : energies_dmft%e_hu_qmc(:) = zero
166 : 106 : energies_dmft%e_dc_tot = zero
167 : 106 : energies_dmft%e_dcdc = zero
168 : 106 : energies_dmft%e_hu_dftu_tot = zero
169 : 106 : energies_dmft%e_hu_mig_tot = zero
170 : 106 : energies_dmft%e_hu_qmc_tot = zero
171 : 106 : energies_dmft%e_hu_tot = zero
172 : 106 : energies_dmft%eband_dft = zero
173 : 106 : energies_dmft%eband_dmft = zero
174 : 106 : energies_dmft%edmft = zero
175 : 106 : energies_dmft%ekin_imp = zero
176 : 106 : energies_dmft%emig_imp = zero
177 : 106 : energies_dmft%emig_loc = zero
178 : 106 : energies_dmft%fband_dft = zero
179 : 106 : energies_dmft%fband_dmft = zero
180 : 106 : energies_dmft%fband_imp = zero
181 : 106 : energies_dmft%fband_weiss = zero
182 : 106 : energies_dmft%fdmft = zero
183 : 106 : energies_dmft%fimp = zero
184 : 106 : energies_dmft%integral = zero
185 : : !energies_dmft%natom = natom
186 : 106 : energies_dmft%sdmft = zero
187 : 106 : energies_dmft%simp = zero
188 : :
189 : 106 : 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 : 106 : 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 : 106 : paw_dmft%e_dc = energies_dmft%e_dc_tot
217 : 106 : paw_dmft%e_hu = energies_dmft%e_hu_tot
218 : 106 : paw_dmft%sdmft = energies_dmft%sdmft
219 : 106 : paw_dmft%simp = energies_dmft%simp
220 : :
221 : 106 : energies_dmft%e_hu => null()
222 [ + - ]: 106 : ABI_SFREE(energies_dmft%e_dc)
223 [ + - ]: 106 : ABI_SFREE(energies_dmft%e_hu_dftu)
224 [ + - ]: 106 : ABI_SFREE(energies_dmft%e_hu_mig)
225 [ + - ]: 106 : ABI_SFREE(energies_dmft%e_hu_qmc)
226 : :
227 : 106 : 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 : 304 : 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 [ + + ]: 304 : if (abs(pawprtvol) >= 3) then
263 [ + + ]: 968 : do iatom=1,paw_dmft%natom
264 : 748 : lpawu = paw_dmft%lpawu(iatom)
265 [ + + ]: 748 : if (lpawu == -1) cycle
266 : 244 : write(tag,'(i4)') iatom
267 [ + - ]: 244 : write(message,'(a,4x,3a)') ch10,"For Correlated Atom ",trim(adjustl(tag)),","
268 : 244 : call wrtout(std_out,message,'COLL')
269 : 244 : write(message,'(26x,a,f12.6)') " E_hu =",energies_dmft%e_hu(iatom)
270 : 244 : call wrtout(std_out,message,'COLL')
271 : 244 : write(message,'(26x,a,f12.6)') " E_hu_mig =",energies_dmft%e_hu_mig(iatom)
272 : 244 : call wrtout(std_out,message,'COLL')
273 : 244 : write(message,'(26x,a,f12.6)') " E_hu_qmc =",energies_dmft%e_hu_qmc(iatom)
274 : 244 : call wrtout(std_out,message,'COLL')
275 : 244 : write(message,'(26x,a,f12.6)') " E_hu_dftu =",energies_dmft%e_hu_dftu(iatom)
276 : 244 : call wrtout(std_out,message,'COLL')
277 : 244 : write(message,'(26x,a,f12.6)') " E_dc =",energies_dmft%e_dc(iatom)
278 : 968 : call wrtout(std_out,message,'COLL')
279 : : end do ! iatom
280 : : end if ! abs(pawprtvol)>=3
281 : 304 : write(message,'(a,5x,2a,5x,a,9(a,5x,a,2x,f15.11),a,5x,a)') ch10, &
282 : 304 : & "-----------------------------------------------",ch10, &
283 : 304 : & "--- Energy in DMFT (in Ha) ",ch10, &
284 : 304 : & "--- E_bandlda (1) (Ha.) = ",energies_dmft%eband_dft,ch10, &
285 : 304 : & "--- E_banddmft(2) (Ha.) = ",energies_dmft%eband_dmft,ch10, &
286 : 304 : & "--- E_hu (3) (Ha.) = ",energies_dmft%e_hu_tot,ch10, &
287 : 304 : & "--- E_hu_mig (4) (Ha.) = ",energies_dmft%e_hu_mig_tot,ch10, &
288 : 304 : & "--- E_hu_qmc (4) (Ha.) = ",energies_dmft%e_hu_qmc_tot,ch10, &
289 : 304 : & "--- E_hu_dftu (5) (Ha.) = ",energies_dmft%e_hu_dftu_tot,ch10, &
290 : 304 : & "--- E_dc (6) (Ha.) = ",energies_dmft%e_dc_tot,ch10, &
291 : 304 : & "--- edmft=( 3-6)(Ha.) = ",energies_dmft%edmft,ch10, &
292 : 304 : & "--- (2-1+3-6)(Ha.) = ",energies_dmft%eband_dmft-energies_dmft%eband_dft+energies_dmft%edmft,ch10, &
293 : 608 : & "-----------------------------------------------"
294 : 304 : call wrtout(std_out,message,'COLL')
295 [ + + ]: 304 : if (idmftloop >= 1) then
296 : 110 : write(message,'(a,i3,1x,f15.11,a)') " (Edmft",idmftloop,energies_dmft%edmft,")"
297 : 110 : call wrtout(ab_out,message,'COLL')
298 : : end if
299 : :
300 : 304 : 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 : 414 : 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 [ + + ]: 414 : if (part == 'both') then
348 [ + + ]: 207 : if (occ_type == " lda") then
349 : 97 : write(message,'(2a)') ch10," == Check: Compute DFT energy terms"
350 : : else
351 : 110 : write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms"
352 : : end if
353 : 207 : call wrtout(std_out,message,'COLL')
354 [ + + ]: 207 : else if (part == 'band') then
355 : 97 : write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms : Band energy terms"
356 : 97 : call wrtout(std_out,message,'COLL')
357 [ - + ]: 110 : 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 [ + - - + ]: 414 : 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 [ + + + + ]: 414 : if (part == 'band' .or. part == 'both') then
370 : 304 : 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 [ + - + + ]: 414 : if (part == 'corr' .or. part == 'both') then
406 : :
407 : : ! == Compute Correlation energy from Migdal formula
408 : : ! -----------------------------------------------------------------------
409 [ + + ]: 207 : if (occ_type /= " lda") then
410 : 110 : 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 : 207 : call compute_dftu_energy(energies_dmft,green,paw_dmft,pawtab(:))
417 [ + + ]: 207 : 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 [ + + - + : 201 : & .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 [ + + - + : 76 : 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 [ + + ]: 76 : if (paw_dmft%dmft_solv == 8 .or. paw_dmft%dmft_solv == 10) then
434 : 8 : write(message,'(2a)') ch10,"Warning, energy is recently computed, not checked"
435 : 8 : call wrtout(std_out,message,'COLL')
436 : : end if
437 : : ! == Compute Correlation energy from QMC correlations.
438 : : ! -----------------------------------------------------------------------
439 : 76 : energies_dmft%e_hu_qmc_tot = zero
440 [ + + ]: 232 : do iatom=1,paw_dmft%natom
441 : 156 : lpawu = paw_dmft%lpawu(iatom)
442 [ + + ]: 156 : if (lpawu == -1) cycle
443 : 92 : energies_dmft%e_hu_qmc(iatom) = green%ecorr_qmc(iatom)
444 : 232 : energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_qmc_tot + energies_dmft%e_hu_qmc(iatom)
445 : : end do ! iatom
446 : 76 : energies_dmft%e_hu => energies_dmft%e_hu_qmc(:)
447 : 76 : 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 : 207 : 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 [ + + ]: 414 : if (part /= 'none') then
456 : 304 : 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 : 414 : 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 : 304 : 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 [ + + ]: 304 : if (occ_type == " lda") then
510 : 97 : write(message,'(2a)') ch10," == Compute DFT Band Energy terms"
511 [ - + ]: 207 : else if (present(fcalc_dft)) then
512 : 0 : write(message,'(2a)') ch10," == Compute DFT Free Energy terms"
513 : : else
514 : 207 : write(message,'(2a)') ch10," == Compute DMFT Band Energy terms"
515 : : end if
516 : 304 : call wrtout(std_out,message,'COLL')
517 : 304 : beta = one / paw_dmft%temp
518 : :
519 : : ! == Compute Band Energy
520 : : ! -----------------------------------------------------------------------
521 [ + + ]: 304 : if (occ_type == " lda") energies_dmft%eband_dft = zero
522 [ - + ]: 304 : 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 : 304 : totch2 = zero
526 : : !totch3 = zero
527 : :
528 : 304 : nkpt = paw_dmft%nkpt
529 : 304 : nspinor = paw_dmft%nspinor
530 : 304 : nsppol = paw_dmft%nsppol
531 : :
532 : 304 : band_index = 0
533 [ + + ]: 734 : do isppol=1,nsppol
534 [ + + ]: 4226 : do ikpt=1,nkpt
535 : 3492 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
536 : 3492 : wtk = paw_dmft%wtk(ikpt)
537 : 3492 : ibc = 0
538 [ + + ]: 104664 : do ib=1,nband_k
539 [ + + - + ]: 101172 : if ((.not. paw_dmft%band_in(ib)) .and. (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7)) cycle
540 [ + - ]: 47588 : if (paw_dmft%band_in(ib)) ibc = ibc + 1
541 : 47588 : eig = paw_dmft%eigen(ib+band_index)
542 [ - + ]: 51080 : 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 [ + + ]: 47588 : if (occ_type == " lda") then ! usual calculation: total non interacting energy
549 : 15292 : 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 [ + - ]: 15292 : if (present(ecalc_dft)) then
555 : 15292 : if (ecalc_dft == 1 .or. ecalc_dft == 3) fermie_used = paw_dmft%fermie_dft
556 [ - + ]: 15292 : if (ecalc_dft == 2 .or. ecalc_dft == 4) fermie_used = paw_dmft%fermie ! only for B3 terms
557 : 15292 : occ = occup_fd(eig,fermie_used,paw_dmft%temp)
558 [ - + ]: 15292 : 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 : 15292 : energies_dmft%eband_dft = energies_dmft%eband_dft + occ*eig*wtk
566 : : end if ! occ_type=" lda"
567 [ + - ]: 47588 : if (paw_dmft%band_in(ib)) then
568 : 47588 : 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 : 47588 : 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 [ + - ]: 47588 : 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 : 3922 : band_index = band_index + nband_k
580 : : end do ! ikpt
581 : : end do ! isppol
582 : :
583 [ - + ]: 304 : 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 [ + + + + ]: 304 : if (occ_type == " lda" .and. nsppol == 1 .and. nspinor == 1) energies_dmft%eband_dft = two * energies_dmft%eband_dft
590 [ + + ]: 304 : 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 [ + - ]: 304 : if (present(ecalc_dft)) then
596 [ - + ]: 304 : 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 [ + + ]: 304 : if (occ_type == " lda") then
601 [ + + ]: 97 : if (abs(energies_dmft%eband_dft-energies_dmft%eband_dmft) > tol5) then
602 : 93 : write(message,'(5x,3a,15x,a,f12.6,a,15x,a,5x,f12.5)') "Warning: ", &
603 : 93 : & "Differences between band energy with Fermi-Dirac occupations",ch10, &
604 : 93 : & "and occupations from DFT Green's function is:",energies_dmft%eband_dft-energies_dmft%eband_dmft,ch10, &
605 : 186 : & "which is larger than",tol5
606 : 93 : call wrtout(std_out,message,'COLL')
607 : : write(message,'(a)') &
608 : 93 : & " Action: increase the number of frequencies, or reduce the number of high energy DMFT bands"
609 : 93 : 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 [ - + ]: 304 : 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 : 304 : 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 : 110 : 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 : 110 : complex(dp), allocatable :: omega_fac(:),trace_moments(:,:),trace(:)
679 : 110 : type(matlu_type), allocatable :: self_nwlo_re(:)
680 : 110 : type(matlu_type), pointer :: matlu_tmp(:) => null()
681 : : character(len=500) :: message
682 : : ! *********************************************************************
683 : :
684 : : ! Only imaginary frequencies here
685 [ + - - + ]: 110 : 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 : 110 : beta = one / paw_dmft%temp
693 : 110 : myproc = paw_dmft%myproc
694 : 110 : natom = paw_dmft%natom
695 : 110 : nmoments = 0
696 : 110 : nspinor = paw_dmft%nspinor
697 : 110 : nsppol = paw_dmft%nsppol
698 : 110 : nwlo = green%nw
699 : 110 : temp = paw_dmft%temp
700 : :
701 [ - + ]: 110 : if (self%has_moments == 1) nmoments = self%nmoments
702 : 110 : iatom_ = 0
703 [ - + ]: 110 : 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 [ - + ]: 110 : 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 [ - + + - : 440 : ABI_MALLOC(trace_moments,(natom,nmoments))
+ - - - +
- - + - +
- + ]
715 [ - + + - : 330 : ABI_MALLOC(trace,(natom))
+ - + - -
+ - + ]
716 : :
717 [ + + ]: 496 : e_hu_migdal(:) = zero
718 [ + + ]: 496 : trace(:) = czero
719 : :
720 [ - + ]: 110 : 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 [ - + + - : 716 : ABI_MALLOC(self_nwlo_re,(natom))
+ - + - -
+ - + +
+ ]
730 [ - + + - : 606 : ABI_MALLOC(matlu_tmp,(natom))
- + - + +
+ ]
731 [ - + ]: 110 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),matlu_tmp(:))
732 : 110 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),self_nwlo_re(:))
733 : 110 : call copy_matlu(self%oper(nwlo)%matlu(:),self_nwlo_re(:),natom,opt_re=1)
734 : : end if ! moments
735 : :
736 [ + + ]: 25995 : do ifreq=1,nwlo
737 : :
738 [ + + ]: 25885 : if (self%distrib%procf(ifreq) /= myproc) cycle
739 : :
740 : 8731 : omega = cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
741 : :
742 [ - + ]: 8731 : if (self%has_moments == 1) then
743 : 0 : matlu_tmp => self%oper(ifreq)%matlu(:)
744 : : else
745 [ - + ]: 8731 : call add_matlu(self%oper(ifreq)%matlu(:),self_nwlo_re(:),matlu_tmp(:),natom,-1)
746 : : end if ! moments
747 : :
748 [ - + ]: 8731 : call trace_prod_matlu(matlu_tmp(:),green%oper(ifreq)%matlu(:),natom,trace(:),iatom=iatom_)
749 : :
750 [ + + ]: 44455 : e_hu_migdal(:) = e_hu_migdal(:) + dble(trace(:))*paw_dmft%wgt_wlo(ifreq)*temp*two
751 : :
752 : : end do ! ifreq
753 : :
754 : 110 : call xmpi_sum(e_hu_migdal(:),paw_dmft%spacecomm,ierr)
755 : :
756 [ + - - - : 330 : ABI_MALLOC(omega_fac,(nmoments))
+ - - + -
+ - + ]
757 : :
758 [ - + ]: 110 : 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 [ - - ]: 110 : e_hu_migdal(:) = e_hu_migdal(:) + dble(trace_moments(:,i)*omega_fac(i))
768 : : end do
769 : :
770 [ + - ]: 110 : if (self%has_moments /= 1) then
771 [ - + ]: 110 : call destroy_matlu(matlu_tmp(:),natom)
772 [ - + ]: 110 : ABI_FREE(matlu_tmp)
773 : : end if
774 : 110 : matlu_tmp => null()
775 : :
776 [ - + ]: 110 : ABI_FREE(omega_fac)
777 : :
778 [ + - ]: 110 : if (self%has_moments == 0) then
779 : 110 : call trace_prod_matlu(self_nwlo_re(:),green%occup%matlu(:),natom,trace(:))
780 : 110 : call destroy_matlu(self_nwlo_re(:),natom)
781 [ + - + + : 496 : ABI_FREE(self_nwlo_re)
- + - + ]
782 [ + + ]: 496 : e_hu_migdal(:) = e_hu_migdal(:) + dble(trace(:))
783 : : end if
784 : :
785 [ - + ]: 110 : ABI_FREE(trace_moments)
786 : :
787 [ + + ]: 496 : e_hu_migdal(:) = half * e_hu_migdal(:) ! E_mig = 1/2 * Tr(Sig*G)
788 [ + + ]: 496 : e_hu_migdal_tot = sum(e_hu_migdal(:))
789 : :
790 [ - + ]: 110 : 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 [ - + - - : 110 : 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 [ + + + - : 216 : 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 : 216 : real(dp), allocatable :: noccmmp(:,:,:,:),nocctot(:)
908 : : ! *********************************************************************
909 : :
910 : : ! - allocations
911 : : ! -----------------------------------------------------------------------
912 : :
913 : 216 : e_dc = zero
914 : 216 : e_dc_for_s = zero
915 : 216 : e_dcdc = zero
916 : 216 : e_ee = zero
917 : 216 : edftumdc = zero
918 : 216 : edftumdcdc = zero
919 : 216 : nsploop = max(paw_dmft%nsppol,paw_dmft%nspinor**2)
920 : 216 : nocc = nsploop
921 : 216 : t2g = (paw_dmft%dmft_t2g == 1)
922 : 216 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
923 : :
924 : 216 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
925 : :
926 : 216 : isppol = 1
927 : 216 : ispinor = 1
928 : 216 : ispinor1 = 1
929 : :
930 [ - + + - : 648 : ABI_MALLOC(nocctot,(nocc))
+ - + - -
+ - + ]
931 : :
932 : : ! - Loop and call to pawuenergy
933 : : ! -----------------------------------------------------------------------
934 [ + + ]: 972 : do iatom=1,paw_dmft%natom
935 : 756 : lpawu = paw_dmft%lpawu(iatom)
936 [ + + ]: 756 : if (lpawu == -1) cycle
937 : 268 : itypat = paw_dmft%typat(iatom)
938 : 268 : lpawu1 = lpawu
939 [ + + ]: 268 : if (t2g .or. x2my2d) lpawu1 = 2
940 : 268 : ndim = 2*lpawu + 1
941 : 268 : ndim1 = 2*lpawu1 + 1
942 : :
943 [ - + + - : 1340 : ABI_MALLOC(noccmmp,(2,ndim1,ndim1,nocc))
- + + - -
+ + - + -
+ - - + -
+ ]
944 [ + + + + : 49945 : noccmmp(:,:,:,:) = zero
+ + + + ]
945 : :
946 : : ! - Setup nocctot and noccmmp
947 : : ! -----------------------------------------------------------------------
948 [ + + ]: 747 : nocctot(:) = zero ! contains nmmp in the n m representation
949 : : ! Begin loop over spin/spinors to initialize noccmmp
950 [ + + ]: 747 : do idijeff=1,nsploop
951 : :
952 [ + + ]: 479 : if (nsploop <= 2) then
953 : : isppol = idijeff
954 [ + - ]: 108 : else if (nsploop == 4) then
955 : 108 : ispinor = spinor_idxs(1,idijeff)
956 : 108 : 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 [ + + ]: 2928 : do im1=1,ndim
963 : 2449 : ims1 = im1
964 : : ! Correct bug in computation of DFT+U energy in the t2g/x2my2d case with TRIQS
965 [ - + ]: 2449 : if (x2my2d .and. dmft_optim) ims1 = 5
966 [ - + ]: 2449 : if (t2g .and. dmft_optim) ims1 = mt2g(im1)
967 [ + + ]: 16511 : do im=1,ndim
968 : 13583 : ims = im
969 [ - + ]: 13583 : if (x2my2d .and. dmft_optim) ims = 5
970 [ - + ]: 13583 : 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 : 13583 : & dble(green%occup%matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol))
974 : : noccmmp(2,ims,ims1,idijeff) = &
975 : 16032 : & 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 [ + + ]: 747 : if (green%has_charge_matlu_solver /= 2) then
982 [ + + ]: 227 : if (nsploop == 4 .and. idijeff <= 2) then
983 [ + + ]: 208 : do im1=1,ndim
984 : : nocctot(1) = nocctot(1) + &
985 : 208 : & dble(green%occup%matlu(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol))
986 : : end do ! im1
987 [ + + ]: 201 : 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 [ + + ]: 252 : if (nsploop == 4 .and. idijeff == 1) then
995 : 14 : nocctot(1) = green%charge_matlu_solver(2,iatom) ! total nb of elec for nspinor=2 is (iatom,2) !!
996 [ + + ]: 238 : 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 : 268 : xe1 = e_dc
1003 : 268 : xe2 = e_ee
1004 : : ! write(std_out,*)" nocctot(1)",nocctot(1),green%charge_matlu_solver(iatom,1)
1005 : 268 : edftumdc = zero
1006 : 268 : edftumdcdc = zero
1007 [ + + ]: 268 : if (present(renorm)) then
1008 : 9 : upawu = one
1009 : 9 : jpawu = renorm(iatom)
1010 : 9 : prt_pawuenergy = 0
1011 : : else
1012 : 259 : upawu = pawtab(itypat)%upawu
1013 : 259 : jpawu = pawtab(itypat)%jpawu
1014 : 259 : 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 : 268 : & e_ee=e_ee,e_dc=e_dc,e_dcdc=e_dcdc,u_dmft=upawu,j_dmft=jpawu,paw_dmft=paw_dmft)
1019 : :
1020 [ - + ]: 268 : 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 : 268 : energies_dmft%e_dc(iatom) = e_dc - xe1
1027 : 268 : energies_dmft%e_hu_dftu(iatom) = e_ee - xe2
1028 : :
1029 [ - + ]: 972 : ABI_FREE(noccmmp)
1030 : : end do ! iatom
1031 : :
1032 [ - + ]: 216 : ABI_FREE(nocctot)
1033 : :
1034 : : ! - gather results
1035 : : ! -----------------------------------------------------------------------
1036 : 216 : energies_dmft%e_dc_tot = e_dc ! this is the only quantity used afterwards.
1037 : 216 : energies_dmft%e_hu_dftu_tot = e_ee
1038 : 216 : energies_dmft%e_dcdc = e_dcdc
1039 [ - + ]: 216 : 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 : 216 : 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 : : !!***
|
