Line data Source code
1 : !!****m* ABINIT/m_entropyDMFT
2 : !! NAME
3 : !! m_entropyDMFT
4 : !!
5 : !! FUNCTION
6 : !! FIXME: add description.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! NOTES
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_entropyDMFT
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_xmpi
30 : use m_dtset
31 :
32 : use m_energies, only : energies_type
33 : use m_splines, only : spline_integrate, spline, splint
34 : use m_pawang, only : pawang_type
35 : use m_pawrad, only : pawrad_type, simp_gen, poisson
36 : use m_pawtab, only : pawtab_type
37 : use m_paw_correlations,only : pawpuxinit
38 : use m_io_tools, only : get_unit
39 : use m_data4entropyDMFT
40 :
41 : implicit none
42 :
43 : private
44 :
45 : public :: entropyDMFT_init
46 : public :: entropyDMFT_destroy
47 : public :: entropyDMFT_nextLambda
48 : public :: entropyDMFT_addIntegrand
49 : public :: entropyDMFT_computeEntropy
50 :
51 : integer, parameter :: E_U0 = 1
52 : integer, parameter :: E_UU = 2
53 : integer, parameter :: E_DIRECT = 1
54 : integer, parameter :: E_DC = 2
55 : integer, parameter :: AC_NOTHING = 0
56 : integer, parameter :: AC_ETOT = 1
57 : character(len=21), parameter :: HDR_NAME = "DATA FOR ETOT DMFT v="
58 :
59 : !!***
60 :
61 : !!****t* m_entropyDMFT/entropyDMFT
62 : !! NAME
63 : !! entropyDMFT
64 : !!
65 : !! FUNCTION
66 : !! This structured datatype contains the necessary data
67 : !!
68 : !! COPYRIGHT
69 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
70 : !! This file is distributed under the terms of the
71 : !! GNU General Public License, see ~abinit/COPYING
72 : !! or http://www.gnu.org/copyleft/gpl.txt .
73 : !!
74 : !! SOURCE
75 :
76 : type, public :: entropyDMFT_t
77 : logical :: isset=.FALSE. ! flag to be sure we are initialize
78 : integer :: spacecomm ! MPI comm
79 : integer :: rank ! rank in the comm
80 : integer :: comm_size ! Number of cpus in the comm
81 : integer :: action ! what to do in gstate
82 : integer :: mylambda ! Current lambda
83 : integer :: natom ! number of atoms
84 : integer :: ncatom ! number of correlated atoms
85 : integer :: ntypat ! number of type of atoms
86 : integer :: nctypat ! number of type of correlated atoms
87 : integer :: nlambda ! number of integration points
88 : integer :: ofile ! unit for file output data
89 : character(len=fnlen) :: filename ! name for the file user readable
90 : character(len=fnlen) :: filedata ! name for the file for restart purposes
91 : character(len=fnlen) :: ofilename ! ofstream prefix
92 : character(len=fnlen) :: ifilename ! ifstream prefix
93 : real(dp) :: temp ! temperature
94 : real(dp) :: entropy0 ! entropy for lambda=0
95 : real(dp) :: energies(2,2) ! internal energy for lambda=0 and 1
96 : integer , allocatable :: index_atom(:) ! index for correlated atoms
97 : integer , allocatable :: index_typat(:) ! index for correlated types
98 : integer , allocatable :: lpawu(:) ! orbital moment to treat (ntypat,2)
99 : integer , allocatable :: typat(:) ! type of each correlated atom
100 : real(dp), allocatable :: U_input(:) ! U from input file (ntypat)
101 : real(dp), allocatable :: J_input(:) ! J from input file (ntypat)
102 : real(dp), allocatable :: lambda(:) ! ilamda
103 : real(dp), allocatable :: docc(:,:,:) ! n_In_j, natom, ilamda
104 : real(dp), allocatable :: e_dc(:,:) ! natom, ilamda
105 : real(dp), allocatable :: uij(:,:) ! uij to compute <uij n_i n_j>, ntypat
106 : end type entropyDMFT_t
107 : !!***
108 :
109 :
110 : contains
111 : !!***
112 :
113 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_init
114 : !! NAME
115 : !! entropyDMFT_init
116 : !!
117 : !! FUNCTION
118 : !! FIXME: add description.
119 : !!
120 : !! COPYRIGHT
121 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
122 : !! This file is distributed under the terms of the
123 : !! GNU General Public License, see ~abinit/COPYING
124 : !! or http://www.gnu.org/copyleft/gpl.txt .
125 : !!
126 : !! INPUTS
127 : !! argin(sizein)=description
128 : !!
129 : !! OUTPUT
130 : !! argout(sizeout)=description
131 : !!
132 : !! SIDE EFFECTS
133 : !!
134 : !! NOTES
135 : !!
136 : !! SOURCE
137 :
138 3 : subroutine entropyDMFT_init(e_t,dt,pawtab,spacecomm,ifilename,ofilename)
139 :
140 : !Arguments ------------------------------------
141 : type(entropyDMFT_t) , intent(inout) :: e_t
142 : type(dataset_type) , intent(in ) :: dt
143 : type(pawtab_type) , intent(in ) :: pawtab(:)
144 : integer , intent(in ) :: spacecomm
145 : character(len=fnlen), intent(in ) :: ifilename
146 : character(len=fnlen), intent(in ) :: ofilename
147 : !Local variables ------------------------------
148 : logical :: doRestart
149 : integer :: natom, ncatom
150 : integer :: iatom, icatom
151 : integer :: nctypat,itypat, ictypat
152 : integer :: ilambda
153 3 : integer, allocatable :: maptypat(:)
154 : character(len=500) :: message
155 :
156 3 : e_t%action = dt%dmft_entropy
157 3 : e_t%mylambda = e_t%action-1 ! should be usefull to start form a given value of lambda
158 : ! -1 added to start from 1 at the first call
159 : ! of nextLambda
160 :
161 3 : if ( e_t%action == AC_NOTHING ) then
162 0 : e_t%isset = .TRUE.
163 : return
164 : endif
165 :
166 3 : e_t%action = AC_ETOT
167 :
168 3 : natom = dt%natom
169 3 : ncatom = dt%natpawu
170 :
171 3 : icatom = 0
172 14 : do iatom=1,natom
173 14 : if ( pawtab(dt%typat(iatom))%lpawu /= -1 ) icatom = icatom + 1
174 : end do
175 :
176 3 : if ( icatom /= ncatom ) ABI_ERROR("Inconsistent number of correlated atoms")
177 :
178 3 : nctypat = 0
179 10 : do itypat=1,dt%ntypat
180 10 : if ( pawtab(itypat)%lpawu /= -1 ) nctypat = nctypat + 1
181 : end do
182 :
183 3 : e_t%natom = natom
184 3 : e_t%ncatom = ncatom
185 3 : e_t%ntypat = dt%ntypat
186 3 : e_t%nctypat = nctypat
187 :
188 3 : if ( dt%dmft_nlambda < 3 ) then
189 0 : write(message,'(2a,i4,2a)') "DMFT must have dmft_nlamda >= 3 to compute entropy", &
190 0 : "whereas its value is dmft_nlambda = ",dt%dmft_nlambda,ch10,&
191 0 : "Action : check you input variable dmft_nlambda"
192 0 : ABI_ERROR(message)
193 : end if
194 :
195 3 : e_t%nlambda = dt%dmft_nlambda
196 :
197 3 : doRestart = .false.
198 3 : if ( e_t%nlambda < (e_t%mylambda+1) ) then
199 0 : ABI_ERROR("Restart calculation of DMFT entropy with a value of dmft_entropy greater than dmft_nlambda")
200 3 : else if ( e_t%mylambda > 0 ) then
201 3 : doRestart = .true.
202 : end if
203 :
204 3 : call entropyDMFT_allocateAll(e_t)
205 :
206 3 : e_t%entropy0 = zero
207 21 : e_t%energies(:,:) = zero
208 :
209 23 : e_t%lambda(:) = (/ (DBLE(ilambda-1)/DBLE(e_t%nlambda-1), ilambda=1,e_t%nlambda) /)
210 :
211 3 : e_t%temp = dt%tsmear
212 :
213 : ! Save each value of U and J for each correlated type
214 3 : ictypat = 1
215 9 : ABI_MALLOC(maptypat,(1:dt%ntypat))
216 10 : do itypat=1,dt%ntypat
217 10 : if ( pawtab(itypat)%lpawu /= -1 ) then
218 3 : e_t%lpawu(ictypat) = pawtab(itypat)%lpawu
219 3 : if(dt%dmft_t2g==1.and.e_t%lpawu(ictypat)==2) then
220 2 : e_t%lpawu(ictypat)=1
221 : end if
222 3 : e_t%U_input(ictypat) = pawtab(itypat)%upawu
223 3 : e_t%J_input(ictypat) = pawtab(itypat)%jpawu
224 3 : e_t%index_typat(ictypat) = itypat
225 3 : maptypat(itypat) = ictypat
226 3 : ictypat = ictypat + 1
227 : end if
228 : end do
229 :
230 : ! Save type local and global of each correlated atom
231 : ! Get the correct value of lpawu for each type
232 3 : icatom = 1
233 6 : do iatom=1,e_t%ncatom
234 6 : if ( pawtab(dt%typat(iatom))%lpawu /= -1 ) then
235 3 : e_t%typat(icatom) = maptypat(dt%typat(iatom))
236 3 : e_t%index_atom(icatom) = iatom
237 3 : icatom = icatom + 1
238 : end if
239 : end do
240 3 : ABI_FREE(maptypat)
241 :
242 231 : write(message,'(a,1x,78a)') ch10,"+",(/ ("-",ilambda=1,76) /), "+"
243 3 : call wrtout(std_out,message,"COLL")
244 3 : call wrtout(ab_out,message,"COLL")
245 3 : write(message,'(1x,a)') "| Calculation of entropy within the DMFT Framework |"
246 3 : call wrtout(std_out,message,"COLL")
247 3 : call wrtout(ab_out,message,"COLL")
248 117 : write(message,'(1x,40a)') "+",(/ ("- ",ilambda=1,38) /), "+"
249 3 : call wrtout(std_out,message,"COLL")
250 13 : do ilambda = 1, e_t%nlambda
251 : write(message,'(1x,a,i4,a11,f6.4,55x,a)') &
252 10 : "|", ilambda, ") lambda = ", e_t%lambda(ilambda), "|"
253 10 : call wrtout(std_out,message,"COLL")
254 23 : do ictypat=1,e_t%nctypat
255 10 : write(message,'(1x,a,6x,a12,i4,4x,a3,5x,2(3x,a4,f6.4,1x,a2),10x,a)') "|", &
256 10 : "- Atom type ", e_t%index_typat(ictypat), "->", &
257 10 : "U = ",e_t%U_input(ictypat)*e_t%lambda(ilambda), "Ha", &
258 20 : "J = ",e_t%J_input(ictypat)*e_t%lambda(ilambda), "Ha","|"
259 20 : call wrtout(std_out,message,"COLL")
260 : end do
261 : end do
262 231 : write(message,'(1x,78a)') "+",(/ ("-",ilambda=1,76) /), "+"
263 3 : call wrtout(std_out,message,"COLL")
264 3 : call wrtout(ab_out,message,"COLL")
265 :
266 : ! Set up MPI
267 3 : e_t%spacecomm = spacecomm
268 3 : e_t%rank = xmpi_comm_rank(spacecomm)
269 3 : e_t%comm_size = xmpi_comm_size(spacecomm)
270 :
271 3 : e_t%ofile = get_unit()
272 3 : e_t%ofilename = ofilename
273 3 : e_t%ifilename = ifilename
274 3 : e_t%filename = TRIM(e_t%ofilename)//"_EntropyDMFT"
275 3 : e_t%filedata = TRIM(e_t%ofilename)//"_data4EntropyDMFT"
276 :
277 3 : if ( doRestart .eqv. .true. ) then
278 : ! If restart fails, then nothing changes and the full calculation is
279 : ! perform. Otherwise, we complete as much a possible the structure.
280 1 : call entropyDMFT_restart(e_t)
281 : end if
282 :
283 : ! Rewrite the files with the previous data
284 3 : call entropyDMFT_dump(e_t)
285 :
286 3 : e_t%isset = .TRUE.
287 :
288 :
289 : end subroutine entropyDMFT_init
290 : !!***
291 :
292 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_allocateAll
293 : !! NAME
294 : !! entropyDMFT_allocateAll
295 : !!
296 : !! FUNCTION
297 : !! FIXME: add description.
298 : !!
299 : !! COPYRIGHT
300 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
301 : !! This file is distributed under the terms of the
302 : !! GNU General Public License, see ~abinit/COPYING
303 : !! or http://www.gnu.org/copyleft/gpl.txt .
304 : !!
305 : !! INPUTS
306 : !! argin(sizein)=description
307 : !!
308 : !! OUTPUT
309 : !! argout(sizeout)=description
310 : !!
311 : !! SIDE EFFECTS
312 : !!
313 : !! NOTES
314 : !!
315 : !! SOURCE
316 :
317 3 : subroutine entropyDMFT_allocateAll(e_t)
318 :
319 : !Arguments ------------------------------------
320 : type(entropyDMFT_t), intent(inout) :: e_t
321 :
322 9 : ABI_MALLOC(e_t%index_atom, (1:e_t%natom))
323 14 : e_t%index_atom = 0
324 9 : ABI_MALLOC(e_t%index_typat,(1:e_t%nctypat))
325 14 : e_t%index_atom = 0
326 9 : ABI_MALLOC(e_t%typat, (1:e_t%ncatom))
327 6 : e_t%typat = 0
328 6 : ABI_MALLOC(e_t%lpawu, (1:e_t%nctypat))
329 6 : e_t%lpawu = 0
330 9 : ABI_MALLOC(e_t%U_input, (1:e_t%nctypat))
331 6 : e_t%U_input = zero
332 6 : ABI_MALLOC(e_t%J_input, (1:e_t%nctypat))
333 6 : e_t%J_input = zero
334 9 : ABI_MALLOC(e_t%lambda, (1:e_t%nlambda))
335 13 : e_t%lambda = 0
336 12 : ABI_MALLOC(e_t%docc, (1:(14*13)/2,1:e_t%ncatom,1:e_t%nlambda))
337 933 : e_t%docc = zero
338 12 : ABI_MALLOC(e_t%e_dc, (1:e_t%ncatom,1:e_t%nlambda))
339 23 : e_t%e_dc = zero
340 9 : ABI_MALLOC(e_t%uij, (1:(14*13)/2,1:e_t%nctypat))
341 279 : e_t%uij = zero
342 3 : end subroutine entropyDMFT_allocateAll
343 : !!***
344 :
345 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_restart
346 : !! NAME
347 : !! entropyDMFT_restart
348 : !!
349 : !! FUNCTION
350 : !! FIXME: add description.
351 : !!
352 : !! COPYRIGHT
353 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
354 : !! This file is distributed under the terms of the
355 : !! GNU General Public License, see ~abinit/COPYING
356 : !! or http://www.gnu.org/copyleft/gpl.txt .
357 : !!
358 : !! INPUTS
359 : !! argin(sizein)=description
360 : !!
361 : !! OUTPUT
362 : !! argout(sizeout)=description
363 : !!
364 : !! SIDE EFFECTS
365 : !!
366 : !! NOTES
367 : !!
368 : !! SOURCE
369 :
370 1 : subroutine entropyDMFT_restart(e_t)
371 :
372 : !Arguments ------------------------------------
373 : type(entropyDMFT_t), intent(inout) :: e_t
374 : !Local variables-------------------------------
375 : logical :: doBcast
376 : character(len=200) :: msg
377 : character(len=21 ) :: hdr
378 : integer :: iostate
379 : integer :: ilambda
380 : integer :: tlambda
381 : integer :: ictypat
382 : integer :: ndim
383 : integer :: icouple
384 : integer :: iflavor1
385 : integer :: iflavor2
386 : integer :: icatom
387 : integer :: iatom
388 : real(dp) :: lambda
389 :
390 1 : write(msg,'(a,1x,2a,1x,a,f6.4,a)') ch10,"EtotDMFT will try to restart calculation from a previous run", ch10, &
391 2 : "and this calculation should start at lambda = ",e_t%lambda(e_t%mylambda+1), ch10
392 1 : call wrtout(std_out,msg,"COLL")
393 1 : call wrtout(ab_out,msg,"COLL")
394 :
395 1 : doBcast = .true.
396 :
397 1 : if ( e_t%rank == 0 ) then
398 1 : tlambda = 1
399 1 : inquire(file=e_t%filedata,iostat=iostate)
400 :
401 1 : if ( iostate /= 0 ) then
402 0 : write(msg,'(5a)') "File ", trim(e_t%filedata), " does not exist or is not accessible", ch10, &
403 0 : "-> No restart performed but full calculation."
404 0 : ABI_WARNING(msg)
405 0 : e_t%mylambda = 0
406 : end if
407 :
408 1 : open(unit=e_t%ofile,file=e_t%filedata,action="read",form="unformatted")
409 1 : read(e_t%ofile,end=42) hdr, iostate
410 :
411 1 : if ( hdr /= HDR_NAME .or. iostate /= 1 ) then
412 0 : write(msg,'(5a)') "File ", trim(e_t%filedata), " does not contain the proper header", ch10, &
413 0 : "-> No restart performed but full calculation."
414 0 : ABI_WARNING(msg)
415 0 : e_t%mylambda = 0
416 : end if
417 :
418 2 : do ilambda = 1, e_t%mylambda
419 2 : tlambda = ilambda
420 2 : read(e_t%ofile,end=42) lambda
421 2 : if ( ABS(lambda - e_t%lambda(ilambda)) >= tol9 ) then
422 1 : write(msg,'(5a,f6.4,a,f6.4)') "File ", trim(e_t%filedata), " is wrong:", ch10, &
423 2 : "Lambda values are different: in file ", lambda, " instead of ", e_t%lambda(ilambda)
424 1 : ABI_WARNING(msg)
425 1 : goto 42
426 : end if
427 :
428 1 : if ( ilambda == 1 ) then
429 1 : read(e_t%ofile,end=42) lambda
430 1 : if ( lambda /= e_t%temp ) then
431 0 : write(msg,'(7a,f6.4)') "File ", trim(e_t%filedata), " is wrong:", ch10, &
432 0 : "Temperature is different than the value of tsmear", ch10, &
433 0 : "-> No restart performed but full calculation."
434 0 : ABI_WARNING(msg)
435 0 : goto 42
436 : end if
437 1 : read(e_t%ofile,end=42) e_t%entropy0
438 1 : read(e_t%ofile,end=42) e_t%energies(E_DC,E_U0)
439 0 : else if ( ilambda == e_t%nlambda ) then ! should never happend ?!
440 0 : read(e_t%ofile,end=42) e_t%energies(E_DC,E_UU)
441 0 : do ictypat = 1, e_t%nctypat
442 0 : ndim = 2*(2*e_t%lpawu(ictypat)+1)
443 0 : icouple = 0
444 0 : do iflavor1 = 1, ndim
445 0 : do iflavor2 = iflavor1+1, ndim
446 0 : icouple = icouple + 1
447 0 : read(e_t%ofile,end=42) e_t%uij(icouple,ictypat)
448 : end do
449 : end do
450 : end do
451 : end if
452 :
453 2 : do icatom = 1, e_t%ncatom
454 1 : iatom = e_t%index_atom(icatom)
455 1 : ndim = 2*(2*e_t%lpawu(e_t%typat(icatom))+1)
456 1 : icouple = 0
457 1 : read(e_t%ofile,end=42) e_t%e_dc(icatom,ilambda)
458 8 : do iflavor1 = 1, ndim
459 22 : do iflavor2 = iflavor1+1, ndim
460 15 : icouple = icouple + 1
461 21 : read(e_t%ofile,end=42) e_t%docc(icouple,icatom,ilambda)
462 : end do
463 : end do
464 : end do
465 : end do
466 0 : close(e_t%ofile)
467 0 : goto 43
468 1 : 42 write(msg,'(5a,f6.4)') "File ", trim(e_t%filedata), " is wrong or incomplete", ch10, &
469 2 : "-> Restart calculation will restart at lambda = ",e_t%lambda(tlambda)
470 1 : ABI_WARNING(msg)
471 1 : close(e_t%ofile)
472 1 : e_t%mylambda = tlambda-1 ! -1 to go to previous lambda
473 : end if
474 : ! MPI BDCAST
475 1 : 43 call xmpi_bcast(e_t%mylambda,0, e_t%spacecomm, ictypat)
476 1 : call xmpi_bcast(e_t%entropy0,0, e_t%spacecomm, ictypat)
477 1 : call xmpi_bcast(e_t%energies,0, e_t%spacecomm, ictypat)
478 : !call xmpi_bcast(e_t%uij,0, e_t%spacecomm, ictypat) ! No need since it
479 : ! restart always perform the last lambda and this is calculation at the end
480 : ! of last lambda
481 1 : call xmpi_bcast(e_t%e_dc,0, e_t%spacecomm, ictypat)
482 1 : call xmpi_bcast(e_t%docc,0, e_t%spacecomm, ictypat)
483 :
484 1 : end subroutine entropyDMFT_restart
485 : !!***
486 :
487 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_dump
488 : !! NAME
489 : !! entropyDMFT_dump
490 : !!
491 : !! FUNCTION
492 : !! FIXME: add description.
493 : !!
494 : !! COPYRIGHT
495 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
496 : !! This file is distributed under the terms of the
497 : !! GNU General Public License, see ~abinit/COPYING
498 : !! or http://www.gnu.org/copyleft/gpl.txt .
499 : !!
500 : !! INPUTS
501 : !! argin(sizein)=description
502 : !!
503 : !! OUTPUT
504 : !! argout(sizeout)=description
505 : !!
506 : !! SIDE EFFECTS
507 : !!
508 : !! NOTES
509 : !!
510 : !! SOURCE
511 :
512 3 : subroutine entropyDMFT_dump(e_t)
513 :
514 : !Arguments ------------------------------------
515 : type(entropyDMFT_t), intent(inout) :: e_t
516 : !Local variables-------------------------------
517 : integer :: ilambda
518 : integer :: ictypat
519 : integer :: ndim
520 : integer :: icouple
521 : integer :: iflavor1
522 : integer :: iflavor2
523 : integer :: icatom
524 : integer :: iatom
525 :
526 3 : if ( e_t%rank /= 0 ) return
527 :
528 : ! Dump _data4EtotDMFT
529 3 : open(unit=e_t%ofile,file=e_t%filedata,action="write",form="unformatted")
530 3 : write(e_t%ofile) HDR_NAME, 1
531 :
532 4 : do ilambda = 1, e_t%mylambda
533 1 : write(e_t%ofile) e_t%lambda(ilambda)
534 :
535 1 : if ( ilambda == 1 ) then
536 1 : write(e_t%ofile) e_t%temp
537 1 : write(e_t%ofile) e_t%entropy0
538 1 : write(e_t%ofile) e_t%energies(E_DC,E_U0)
539 0 : else if ( ilambda == e_t%nlambda ) then ! should never happend ?!
540 0 : write(e_t%ofile) e_t%energies(E_DC,E_UU)
541 0 : do ictypat = 1, e_t%nctypat
542 0 : ndim = 2*(2*e_t%lpawu(ictypat)+1)
543 0 : icouple = 0
544 0 : do iflavor1 = 1, ndim
545 0 : do iflavor2 = iflavor1+1, ndim
546 0 : icouple = icouple + 1
547 0 : write(e_t%ofile) e_t%uij(icouple,ictypat)
548 : end do
549 : end do
550 : end do
551 : end if
552 :
553 5 : do icatom = 1, e_t%ncatom
554 1 : iatom = e_t%index_atom(icatom)
555 1 : ndim = 2*(2*e_t%lpawu(e_t%typat(icatom))+1)
556 1 : icouple = 0
557 1 : write(e_t%ofile) e_t%e_dc(icatom,ilambda)
558 8 : do iflavor1 = 1, ndim
559 22 : do iflavor2 = iflavor1+1, ndim
560 15 : icouple = icouple + 1
561 21 : write(e_t%ofile) e_t%docc(icouple,icatom,ilambda)
562 : end do
563 : end do
564 : end do
565 : end do
566 3 : close(e_t%ofile)
567 :
568 : ! Dump _EtotDMFT
569 3 : open(unit=e_t%ofile,file=e_t%filename)
570 3 : write(e_t%ofile,'(2a)') "# Data for entropy calculation in DMFT",ch10
571 4 : do ilambda = 1, e_t%mylambda
572 4 : if ( ilambda == 1 ) then
573 1 : write(e_t%ofile,'(a)') "# Temperature [Ha]:"
574 1 : write(e_t%ofile,'(es22.14)') e_t%temp
575 1 : write(e_t%ofile,'(a)') "# Entropy for lambda=0 [kb]:"
576 1 : write(e_t%ofile,'(es22.14)') e_t%entropy0
577 1 : write(e_t%ofile,'(a)') "# Internal energy for lambda=0 [Ha]:"
578 1 : write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_U0)
579 0 : else if ( e_t%mylambda == e_t%nlambda ) then
580 0 : write(e_t%ofile,'(a)') "# Internal energy for lambda=1 [Ha]:"
581 0 : write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_UU)
582 0 : do ictypat = 1, e_t%nctypat
583 0 : ndim = 2*(2*e_t%lpawu(ictypat)+1)
584 0 : icouple = 0
585 0 : write(e_t%ofile,'(a,f7.5,1x,a,i4)') "# Interaction Matrix normalized by U=",e_t%U_input(ictypat) , &
586 0 : "[Ha] for atom type", e_t%index_typat(ictypat)
587 0 : do iflavor1 = 1, ndim
588 0 : write(e_t%ofile,'(14(a21,2x))',advance="no") (/ ( "...", iflavor2=1,iflavor1 ) /)
589 0 : do iflavor2 = iflavor1+1, ndim
590 0 : icouple = icouple + 1
591 0 : write(e_t%ofile,'(14(es21.14,2x))',advance="no") e_t%uij(icouple,ictypat)
592 : end do
593 0 : write(e_t%ofile,*)
594 : end do
595 : end do
596 : end if
597 : end do
598 3 : close(e_t%ofile)
599 : end subroutine entropyDMFT_dump
600 : !!***
601 :
602 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_nextLambda
603 : !! NAME
604 : !! entropyDMFT_nextLambda
605 : !!
606 : !! FUNCTION
607 : !! FIXME: add description.
608 : !!
609 : !! COPYRIGHT
610 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
611 : !! This file is distributed under the terms of the
612 : !! GNU General Public License, see ~abinit/COPYING
613 : !! or http://www.gnu.org/copyleft/gpl.txt .
614 : !!
615 : !! INPUTS
616 : !! argin(sizein)=description
617 : !!
618 : !! OUTPUT
619 : !! argout(sizeout)=description
620 : !!
621 : !! SIDE EFFECTS
622 : !!
623 : !! NOTES
624 : !!
625 : !! SOURCE
626 :
627 12 : function entropyDMFT_nextLambda(e_t,dt,pawtab,pawang,pawrad) result(nextstep)
628 :
629 : !Arguments ------------------------------------
630 : type(entropyDMFT_t) , intent(inout) :: e_t
631 : type(dataset_type) , intent(in) :: dt
632 : type(pawtab_type) , intent(inout) :: pawtab(:)
633 : type(pawang_type) , intent(in ) :: pawang
634 : type(pawrad_type) , intent(inout) :: pawrad(:)
635 : !Local variables ------------------------------
636 : logical :: nextstep
637 : integer :: itypat
638 : integer :: mylambda
639 : logical :: is_dfpt=.false.
640 12 : real(dp),allocatable :: upawu(:),jpawu(:)
641 : character(len=100) :: message
642 :
643 12 : if ( e_t%isset .eqv. .FALSE. ) &
644 0 : ABI_ERROR("entropyDMFT is not initialized")
645 :
646 : ! go to next lambda
647 12 : mylambda = e_t%mylambda + 1
648 : !if ( present(ilambda) ) mylambda = ilambda
649 12 : e_t%mylambda = mylambda
650 :
651 12 : if ( e_t%action == AC_NOTHING .and. mylambda >= 1 ) then
652 : nextstep = .FALSE.
653 : ! We do nothing and return false, one scfvc has already performed
654 12 : else if ( e_t%action == AC_NOTHING .and. mylambda < 1 ) then
655 : nextstep = .TRUE.
656 : ! We do nothing but return true to perform scfcv a least once
657 12 : else if ( e_t%action == AC_ETOT .and. mylambda <= e_t%nlambda ) then ! we iterate over lambda
658 9 : nextstep = .TRUE.
659 27 : ABI_MALLOC(upawu,(dt%ntypat))
660 18 : ABI_MALLOC(jpawu,(dt%ntypat))
661 18 : do itypat = 1, e_t%nctypat
662 9 : upawu(e_t%index_typat(itypat)) = e_t%lambda(mylambda) * e_t%U_input(itypat)
663 18 : jpawu(e_t%index_typat(itypat)) = e_t%lambda(mylambda) * e_t%J_input(itypat)
664 : end do
665 : else ! we did all lambda values
666 : nextstep = .FALSE.
667 : endif
668 :
669 12 : if ( e_t%action == AC_ETOT .and. nextstep .eqv. .true. ) then
670 693 : write(message,'(a,1x,78a)') ch10,"+",(/ ("-",itypat=1,76) /), "+"
671 9 : call wrtout(std_out,message,"COLL")
672 9 : call wrtout(ab_out,message,"COLL")
673 : write(message,'(1x,a,i4,a11,f6.4,55x,a)') &
674 9 : "|", mylambda, ") lambda = ", e_t%lambda(mylambda), "|"
675 9 : call wrtout(std_out,message,"COLL")
676 9 : call wrtout(ab_out,message,"COLL")
677 18 : do itypat=1,e_t%nctypat
678 9 : write(message,'(1x,a,6x,a12,i4,4x,a3,5x,2(3x,a4,f6.4,1x,a2),10x,a)') "|",&
679 9 : "- Atom type ", e_t%index_typat(itypat), "->", &
680 9 : "U = ",e_t%U_input(itypat)*e_t%lambda(mylambda), "Ha", &
681 18 : "J = ",e_t%J_input(itypat)*e_t%lambda(mylambda), "Ha","|"
682 9 : call wrtout(std_out,message,"COLL")
683 18 : call wrtout(ab_out,message,"COLL")
684 : end do
685 693 : write(message,'(1x,78a)') "+",(/ ("-",itypat=1,76) /), "+"
686 9 : call wrtout(std_out,message,"COLL")
687 9 : call wrtout(ab_out,message,"COLL")
688 : call pawpuxinit(dt%dmatpuopt,dt%exchmix,dt%f4of2_sla,dt%f6of2_sla,&
689 : & is_dfpt,jpawu,dt%lexexch,dt%lpawu,dt%nspinor,dt%ntypat,dt%optdcmagpawu,pawang,dt%pawprtvol,&
690 : & pawrad,pawtab,upawu,dt%usedmft,dt%useexexch,dt%usepawu,dmft_orbital=dt%dmft_orbital, &
691 : & dmft_dc=dt%dmft_dc,dmft_orbital_filepath=dt%dmft_orbital_filepath,dmft_yukawa_param=dt%dmft_yukawa_param, &
692 9 : & dmft_yukawa_lambda=dt%dmft_yukawa_lambda,dmft_yukawa_epsilon=dt%dmft_yukawa_epsilon)
693 9 : ABI_FREE(upawu)
694 9 : ABI_FREE(jpawu)
695 : end if
696 :
697 12 : end function entropyDMFT_nextLambda
698 : !!***
699 :
700 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_addIntegrand
701 : !! NAME
702 : !! entropyDMFT_addIntegrand
703 : !!
704 : !! FUNCTION
705 : !! FIXME: add description.
706 : !!
707 : !! COPYRIGHT
708 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
709 : !! This file is distributed under the terms of the
710 : !! GNU General Public License, see ~abinit/COPYING
711 : !! or http://www.gnu.org/copyleft/gpl.txt .
712 : !!
713 : !! INPUTS
714 : !! argin(sizein)=description
715 : !!
716 : !! OUTPUT
717 : !! argout(sizeout)=description
718 : !!
719 : !! SIDE EFFECTS
720 : !!
721 : !! NOTES
722 : !!
723 : !! SOURCE
724 :
725 9 : subroutine entropyDMFT_addIntegrand(e_t,dt,energies,data4etot)
726 :
727 : !Arguments ------------------------------------
728 : type(entropyDMFT_t) , intent(inout) :: e_t
729 : type(dataset_type) , intent(in ) :: dt
730 : type(energies_type), intent(in ) :: energies
731 : type(data4entropyDMFT_t) , intent(in ) :: data4etot
732 : !Local variables ------------------------------
733 : integer :: optdc
734 : integer :: ictypat
735 : integer :: iatom, icatom
736 : integer :: iflavor1, iflavor2, ndim
737 : integer :: icouple
738 :
739 9 : if ( e_t%action == AC_NOTHING ) return
740 :
741 : ! Write lambda for restart
742 9 : if ( e_t%rank == 0 ) then
743 9 : open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
744 9 : write(e_t%ofile) e_t%lambda(e_t%mylambda)
745 9 : close(e_t%ofile)
746 : end if
747 :
748 9 : if ( e_t%mylambda == 1 ) then
749 : ! Save entropy and internal energy for U=0
750 2 : e_t%entropy0 = energies%entropy_ks
751 : ! 1 is for usepaw that is 1 in DMFT, optdc is to know if the DC scheme is
752 : ! calculated.
753 2 : call energies%eval_eint(dt,1,optdc,e_t%energies(E_DIRECT,E_U0),e_t%energies(E_DC,E_U0))
754 2 : if ( e_t%rank == 0 ) then
755 2 : open(unit=e_t%ofile,file=e_t%filename,position="append")
756 2 : write(e_t%ofile,'(a)') "# Temperature [Ha]:"
757 2 : write(e_t%ofile,'(es22.14)') e_t%temp
758 2 : write(e_t%ofile,'(a)') "# Entropy for lambda=0 [kb]:"
759 2 : write(e_t%ofile,'(es22.14)') e_t%entropy0
760 2 : write(e_t%ofile,'(a)') "# Internal energy for lambda=0 [Ha]:"
761 2 : write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_U0)
762 2 : close(e_t%ofile)
763 2 : open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
764 2 : write(e_t%ofile) e_t%temp
765 2 : write(e_t%ofile) e_t%entropy0
766 2 : write(e_t%ofile) e_t%energies(E_DC,E_U0)
767 2 : close(e_t%ofile)
768 : end if
769 7 : else if ( e_t%mylambda == e_t%nlambda ) then
770 : ! Save internal energy for U=Umax
771 3 : call energies%eval_eint(dt,1,optdc,e_t%energies(E_DIRECT,E_UU),e_t%energies(E_DC,E_UU))
772 3 : if ( e_t%rank == 0 ) then
773 3 : open(unit=e_t%ofile,file=e_t%filename,position="append")
774 3 : write(e_t%ofile,'(a)') "# Internal energy for lambda=1 [Ha]:"
775 3 : write(e_t%ofile,'(es22.14)') e_t%energies(E_DC,E_UU)
776 : end if
777 : ! Compute U_max, and uij
778 6 : do ictypat = 1, e_t%nctypat
779 3 : ndim = 2*(2*e_t%lpawu(ictypat)+1)
780 3 : icouple = 0
781 : !write(*,*) matU
782 3 : if ( e_t%rank == 0 ) &
783 3 : write(e_t%ofile,'(a,f7.5,1x,a,i4)') "# Interaction Matrix normalized by U=",e_t%U_input(ictypat) , &
784 6 : "[Ha] for atom type", e_t%index_typat(ictypat)
785 32 : do iflavor1 = 1, ndim
786 26 : if ( e_t%rank == 0 ) &
787 320 : write(e_t%ofile,'(14(a21,2x))',advance="no") (/ ( "...", iflavor2=1,iflavor1 ) /)
788 147 : do iflavor2 = iflavor1+1, ndim
789 121 : icouple = icouple + 1
790 : e_t%uij(icouple,ictypat) = data4etot%hu_dens(iflavor1,iflavor2,e_t%index_typat(ictypat)) &
791 121 : /e_t%U_input(ictypat) ! to modify, this is just the idea
792 121 : if ( e_t%rank == 0 ) &
793 147 : write(e_t%ofile,'(14(es21.14,2x))',advance="no") e_t%uij(icouple,ictypat)
794 : end do
795 29 : if ( e_t%rank == 0 ) write(e_t%ofile,*)
796 : end do
797 : end do
798 :
799 3 : if ( e_t%rank == 0 ) then
800 3 : close(e_t%ofile)
801 3 : open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
802 3 : write(e_t%ofile) e_t%energies(E_DC,E_UU)
803 6 : do ictypat = 1, e_t%nctypat
804 3 : ndim = 2*(2*e_t%lpawu(ictypat)+1)
805 3 : icouple = 0
806 32 : do iflavor1 = 1, ndim
807 150 : do iflavor2 = iflavor1+1, ndim
808 121 : icouple = icouple + 1
809 147 : write(e_t%ofile) e_t%uij(icouple,ictypat)
810 : end do
811 : end do
812 : end do
813 3 : close(e_t%ofile)
814 : end if
815 : endif
816 :
817 : ! For all lambda
818 : ! Save Docc, Nup and Ndwn
819 9 : if ( e_t%rank == 0 ) then
820 9 : open(unit=e_t%ofile,file=e_t%filedata,position="append",form="unformatted")
821 : end if
822 18 : do icatom = 1, e_t%ncatom
823 9 : iatom = e_t%index_atom(icatom)
824 9 : ndim = 2*(2*e_t%lpawu(e_t%typat(icatom))+1)
825 9 : icouple = 0
826 9 : e_t%e_dc(icatom,e_t%mylambda) = data4etot%e_dc(iatom)
827 9 : if ( e_t%rank == 0 ) then
828 9 : write(e_t%ofile) e_t%e_dc(icatom,e_t%mylambda)
829 : end if
830 96 : do iflavor1 = 1, ndim
831 450 : do iflavor2 = iflavor1+1, ndim
832 363 : icouple = icouple + 1
833 363 : e_t%docc(icouple,icatom,e_t%mylambda) = data4etot%Docc(iflavor1,iflavor2,iatom)
834 441 : if ( e_t%rank == 0 ) then
835 363 : write(e_t%ofile) e_t%docc(icouple,icatom,e_t%mylambda)
836 : end if
837 : end do
838 : end do
839 : end do
840 :
841 9 : if ( e_t%rank == 0 ) then
842 9 : close(e_t%ofile)
843 : end if
844 :
845 : !call entropyDMFT_computeIntegrand(e_t)
846 :
847 : end subroutine entropyDMFT_addIntegrand
848 : !!***
849 :
850 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_computeEntropy
851 : !! NAME
852 : !! entropyDMFT_computeEntropy
853 : !!
854 : !! FUNCTION
855 : !! FIXME: add description.
856 : !!
857 : !! COPYRIGHT
858 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
859 : !! This file is distributed under the terms of the
860 : !! GNU General Public License, see ~abinit/COPYING
861 : !! or http://www.gnu.org/copyleft/gpl.txt .
862 : !!
863 : !! INPUTS
864 : !! argin(sizein)=description
865 : !!
866 : !! OUTPUT
867 : !! argout(sizeout)=description
868 : !!
869 : !! SIDE EFFECTS
870 : !!
871 : !! NOTES
872 : !!
873 : !! SOURCE
874 :
875 3 : subroutine entropyDMFT_computeEntropy(e_t,entropy)
876 :
877 : !Arguments ------------------------------------
878 : type(entropyDMFT_t), intent(inout) :: e_t
879 : real(dp) , intent(inout) :: entropy !inout in case we do nothing, avoid to change entropy in NaN
880 : !Local variables ------------------------------
881 : integer :: ilambda
882 : integer :: nlambda
883 : integer :: icatom
884 : integer :: ncatom
885 : integer :: icouple
886 : integer :: ndim
887 3 : real(dp), allocatable :: integrand(:,:)
888 : real(dp) :: docc
889 : real(dp) :: integral
890 : real(dp) :: entropyDirect
891 : character(len=500) :: string
892 : character(len=50) :: i2str
893 :
894 3 : if ( e_t%action == AC_NOTHING ) return
895 :
896 3 : nlambda = e_t%nlambda
897 3 : ncatom = e_t%ncatom
898 :
899 12 : ABI_MALLOC(integrand,(1:nlambda,1:ncatom))
900 16 : integrand(1:nlambda,1:ncatom) = zero
901 :
902 3 : if ( e_t%rank == 0 ) then
903 3 : open(unit=e_t%ofile,file=e_t%filename,position="append")
904 3 : write(e_t%ofile,'(2a)') ch10,"#Decomposition for each lambda per atom"
905 : end if
906 13 : do ilambda = 1, nlambda
907 10 : if ( e_t%rank == 0 ) then
908 10 : write(e_t%ofile,'(a,i2,a)') "# ------| lambda ",ilambda," |------ #"
909 10 : write(e_t%ofile,'(a5,a23,a23)') "#Atom","DC","<uij*n_i*n_j>"
910 : end if
911 23 : do icatom = 1, ncatom
912 10 : integrand(ilambda,icatom) = -e_t%e_dc(icatom,ilambda)
913 10 : ndim = 2*(2*e_t%lpawu(e_t%typat(e_t%index_atom(icatom)))+1) ! nflavors
914 10 : ndim = ndim*(ndim-1)/2 ! number of couples
915 10 : docc = zero
916 388 : do icouple = 1, ndim
917 388 : docc = docc + e_t%uij(icouple,e_t%typat(icatom)) * e_t%docc(icouple,icatom,ilambda)
918 : end do
919 10 : integrand(ilambda,icatom) = integrand(ilambda,icatom) + docc
920 10 : if ( e_t%rank == 0 ) &
921 20 : write(e_t%ofile,'(i5,2es23.14)') e_t%index_atom(icatom),e_t%e_dc(icatom,ilambda),docc
922 : end do
923 : end do
924 :
925 : ! print for test purpose
926 3 : if ( e_t%rank == 0 ) then
927 3 : write(i2str,'(i6)') e_t%ncatom
928 3 : string='(a,a7,'//TRIM(ADJUSTL(i2str))//'(12x,"Atom",i6))'
929 : !open(unit=e_t%ofile,file=e_t%filename,position="append")
930 9 : write(e_t%ofile,string) ch10,"#lambda", (/ (e_t%index_atom(icatom),icatom=1,ncatom) /)
931 3 : string='(1x,f6.4,'//TRIM(ADJUSTL(i2str))//'(1x,e21.14))'
932 13 : do ilambda = 1, nlambda
933 33 : write(e_t%ofile,string) e_t%lambda(ilambda), (/ (integrand(ilambda,icatom),icatom=1,ncatom) /)
934 : end do
935 3 : close(e_t%ofile)
936 : end if
937 :
938 : ! Integrate the integrand for all correlated atoms
939 3 : call entropyDMFT_integrate(e_t,integrand,integral)
940 3 : ABI_FREE(integrand)
941 :
942 231 : write(string,'(a,1x,78a)') ch10,"+",(/ ("-",ilambda=1,76) /), "+"
943 3 : call wrtout(std_out,string,"COLL")
944 3 : call wrtout(ab_out,string,"COLL")
945 3 : write(string,'(1x,a)') "| Calculation of entropy within the DMFT Framework |"
946 3 : call wrtout(std_out,string,"COLL")
947 3 : call wrtout(ab_out,string,"COLL")
948 117 : write(string,'(1x,40a)') "+",(/ ("- ",ilambda=1,38) /), "+"
949 3 : call wrtout(std_out,string,"COLL")
950 :
951 3 : entropyDirect = e_t%entropy0 + ( e_t%energies(E_DIRECT,E_UU) - e_t%energies(E_DIRECT,E_U0) - integral ) /e_t%temp
952 3 : entropy = e_t%entropy0 + ( e_t%energies(E_DC,E_UU) - e_t%energies(E_DC,E_U0) - integral ) /e_t%temp
953 :
954 3 : write(i2str,'(a)') '(1x,a,19x,a11,es21.14,1x,a4,20x,a)'
955 3 : write(string,i2str) "|","Integral = ", integral, "[Ha]", "|"
956 3 : call wrtout(std_out,string,"COLL")
957 3 : write(string,i2str) "|","E(0) = ", e_t%energies(E_DC,E_U0), "[Ha]", "|"
958 3 : call wrtout(std_out,string,"COLL")
959 3 : write(string,i2str) "|","E(U) = ", e_t%energies(E_DC,E_UU), "[Ha]", "|"
960 3 : call wrtout(std_out,string,"COLL")
961 3 : write(string,i2str) "|","S(0) = ", e_t%entropy0, "[kb]", "|"
962 3 : call wrtout(std_out,string,"COLL")
963 3 : write(string,i2str) "|","S(U) = ", entropy, "[kb]", "|"
964 3 : call wrtout(std_out,string,"COLL")
965 117 : write(string,'(1x,40a)') "+",(/ ("- ",ilambda=1,38) /), "+"
966 3 : call wrtout(std_out,string,"COLL")
967 :
968 3 : write(string,'(1x,a,16x,a22,es21.14,1x,a4,12x,a)') "|","-kT*entropy is set to ", -e_t%temp*entropy, "[Ha]", "|"
969 3 : call wrtout(std_out,string,"COLL")
970 3 : call wrtout(ab_out,string,"COLL")
971 :
972 231 : write(string,'(1x,78a)') "+",(/ ("-",ilambda=1,76) /), "+"
973 3 : call wrtout(std_out,string,"COLL")
974 3 : call wrtout(ab_out,string,"COLL")
975 :
976 3 : if ( entropy < zero ) then
977 1 : write(string,'(3a)') "Entropy is negative !!!!",ch10,&
978 2 : "It does not make any sense"
979 1 : ABI_WARNING(string)
980 : end if
981 :
982 3 : if ( abs(entropy-entropyDirect) >= tol3 ) then
983 3 : write(string,'(1x,a,1x,f8.3,1x,a,1x,f8.3,2a)') "Difference between Direct and DC entropies is", abs(entropy-entropyDirect), &
984 6 : "which is greater than", tol3,ch10,"Action : converge better the DMFT and/or DFT loops"
985 3 : ABI_WARNING(string)
986 : end if
987 :
988 :
989 :
990 6 : end subroutine entropyDMFT_computeEntropy
991 : !!***
992 :
993 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_integrate
994 : !! NAME
995 : !! entropyDMFT_integrate
996 : !!
997 : !! FUNCTION
998 : !! FIXME: add description.
999 : !!
1000 : !! COPYRIGHT
1001 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
1002 : !! This file is distributed under the terms of the
1003 : !! GNU General Public License, see ~abinit/COPYING
1004 : !! or http://www.gnu.org/copyleft/gpl.txt .
1005 : !!
1006 : !! INPUTS
1007 : !! argin(sizein)=description
1008 : !!
1009 : !! OUTPUT
1010 : !! argout(sizeout)=description
1011 : !!
1012 : !! SIDE EFFECTS
1013 : !!
1014 : !! NOTES
1015 : !!
1016 : !! SOURCE
1017 :
1018 3 : subroutine entropyDMFT_integrate(e_t,integrand,integral)
1019 :
1020 : !Arguments ------------------------------------
1021 : type(entropyDMFT_t), intent(inout) :: e_t
1022 : real(dp) , intent(in ) :: integrand(:,:)
1023 : real(dp) , intent( out) :: integral
1024 : !Local variables ------------------------------
1025 3 : real(dp), allocatable :: ypp(:)
1026 3 : real(dp), allocatable :: fitx(:)
1027 3 : real(dp), allocatable :: fity(:)
1028 : real(dp) :: integral1
1029 : !real(dp) :: integral2
1030 : !real(dp) :: dx
1031 : integer :: unit
1032 : integer :: Nfit
1033 : integer :: icatom
1034 : integer :: i
1035 :
1036 3 : Nfit = 1000
1037 :
1038 9 : ABI_MALLOC(ypp,(1:e_t%nlambda))
1039 3 : ABI_MALLOC(fitx,(1:Nfit))
1040 3 : ABI_MALLOC(fity,(1:Nfit))
1041 :
1042 3 : integral = zero
1043 :
1044 3 : unit = get_unit()
1045 3 : if ( e_t%rank == 0 ) open(unit=unit,file="fit.log")
1046 6 : do icatom = 1, e_t%ncatom
1047 3 : integral1 = zero
1048 : !integral2 = zero
1049 : !dx = e_t%U_input(icatom)/dble(e_t%nlambda-1)
1050 :
1051 : call spline(e_t%lambda,integrand(:,icatom),e_t%nlambda,&
1052 : (integrand(2,icatom)-integrand(1,icatom))/e_t%lambda(2),& ! first derivative left
1053 : (integrand(e_t%nlambda,icatom)-integrand(e_t%nlambda-1,icatom))/e_t%lambda(2),& ! first derivative right
1054 3 : ypp)
1055 6006 : fitx(1:Nfit)=(/ (dble(i-1)/dble(Nfit-1),i=1,Nfit) /)
1056 3 : call splint(e_t%nlambda,e_t%lambda,integrand(:,icatom),ypp,Nfit,fitx,fity)
1057 3003 : integral1 = (sum(fity(1:Nfit))-(fity(1)+fity(Nfit))*half)*e_t%U_input(e_t%typat(icatom))/dble(Nfit-1)
1058 3 : if ( e_t%rank == 0 ) then
1059 3003 : do i=1,Nfit
1060 3003 : write(unit,'(2ES21.14)') fitx(i),fity(i)
1061 : end do
1062 : end if
1063 : !Unfortunately I don't trust this function and I don't have time to understand it.
1064 : !call spline_integrate(integral2,Nfit,dx,integrand(:,icatom))
1065 : !write(*,*) integral1, integral2
1066 :
1067 : !if ( abs(integral2-integral1) >= tol6 ) then
1068 : ! write(msg,'(1x,a,1x,f8.6,1x,a,1x,f8.6,1x,a,i4)') "Difference between two different ways of integration is", abs(integral2-integral1), &
1069 : ! "which is greater than", tol6, "for correlated atom",e_t%index_atom(icatom)
1070 : ! ABI_WARNING(msg)
1071 : !end if
1072 :
1073 6 : integral = integral+integral1
1074 : end do
1075 3 : if ( e_t%rank == 0 ) close(unit)
1076 :
1077 3 : ABI_FREE(ypp)
1078 3 : ABI_FREE(fitx)
1079 3 : ABI_FREE(fity)
1080 :
1081 3 : end subroutine entropyDMFT_integrate
1082 : !!***
1083 :
1084 : !!****f* ABINIT/m_entropyDMFT/entropyDMFT_destroy
1085 : !! NAME
1086 : !! entropyDMFT_destroy
1087 : !!
1088 : !! FUNCTION
1089 : !! FIXME: add description.
1090 : !!
1091 : !! COPYRIGHT
1092 : !! Copyright (C) 2014-2026 ABINIT group (J. Bieder)
1093 : !! This file is distributed under the terms of the
1094 : !! GNU General Public License, see ~abinit/COPYING
1095 : !! or http://www.gnu.org/copyleft/gpl.txt .
1096 : !!
1097 : !! INPUTS
1098 : !! argin(sizein)=description
1099 : !!
1100 : !! OUTPUT
1101 : !! argout(sizeout)=description
1102 : !!
1103 : !! SIDE EFFECTS
1104 : !!
1105 : !! NOTES
1106 : !!
1107 : !! SOURCE
1108 :
1109 3 : subroutine entropyDMFT_destroy(e_t)
1110 :
1111 : !Arguments ------------------------------------
1112 : type(entropyDMFT_t), intent(inout) :: e_t
1113 :
1114 3 : if ( allocated(e_t%index_atom) ) then
1115 3 : ABI_FREE(e_t%index_atom)
1116 : endif
1117 3 : if ( allocated(e_t%index_typat)) then
1118 3 : ABI_FREE(e_t%index_typat)
1119 : endif
1120 3 : if ( allocated(e_t%typat ) ) then
1121 3 : ABI_FREE(e_t%typat )
1122 : endif
1123 3 : if ( allocated(e_t%lpawu ) ) then
1124 3 : ABI_FREE(e_t%lpawu)
1125 : endif
1126 3 : if ( allocated(e_t%U_input ) ) then
1127 3 : ABI_FREE(e_t%U_input)
1128 : endif
1129 3 : if ( allocated(e_t%J_input ) ) then
1130 3 : ABI_FREE(e_t%J_input)
1131 : endif
1132 3 : if ( allocated(e_t%lambda ) ) then
1133 3 : ABI_FREE(e_t%lambda)
1134 : endif
1135 3 : if ( allocated(e_t%docc ) ) then
1136 3 : ABI_FREE(e_t%docc)
1137 : endif
1138 3 : if ( allocated(e_t%e_dc ) ) then
1139 3 : ABI_FREE(e_t%e_dc)
1140 : endif
1141 3 : if ( allocated(e_t%uij ) ) then
1142 3 : ABI_FREE(e_t%uij)
1143 : endif
1144 3 : e_t%isset = .FALSE.
1145 3 : end subroutine entropyDMFT_destroy
1146 : !!***
1147 0 : end module m_entropyDMFT
1148 : !!***
|
