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