Line data Source code
1 : !!****m* ABINIT/m_hubbard_one
2 : !! NAME
3 : !! m_hubbard_one
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! Solve Anderson model with the density/density Hubbard one approximation
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! INPUTS
16 : !!
17 : !! OUTPUT
18 : !!
19 : !! SOURCE
20 :
21 : #if defined HAVE_CONFIG_H
22 : #include "config.h"
23 : #endif
24 :
25 :
26 : #include "abi_common.h"
27 :
28 : ! nvtx related macro definition
29 : #include "nvtx_macros.h"
30 :
31 : MODULE m_hubbard_one
32 :
33 :
34 : use defs_basis
35 :
36 : #ifdef HAVE_GPU_MARKERS
37 : use m_nvtx_data
38 : #endif
39 :
40 : implicit none
41 :
42 : private
43 :
44 : public :: hubbard_one
45 : !!***
46 :
47 : contains
48 :
49 : !!****f* m_hubbard_one/hubbard_one
50 : !! NAME
51 : !! hubbard_one
52 : !!
53 : !! FUNCTION
54 : !! Solve the hubbard one approximation
55 : !!
56 : !! COPYRIGHT
57 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
58 : !! This file is distributed under the terms of the
59 : !! GNU General Public License, see ~abinit/COPYING
60 : !! or http://www.gnu.org/copyleft/gpl.txt .
61 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
62 : !!
63 : !! INPUTS
64 : !! cryst_struc
65 : !! istep = step of iteration for DFT.
66 : !! dft_occup
67 : !! mpi_enreg=information about MPI parallelization
68 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
69 : !!
70 : !! OUTPUT
71 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
72 : !!
73 : !! NOTES
74 : !!
75 : !! SOURCE
76 :
77 138 : subroutine hubbard_one(cryst_struc,green,hu,paw_dmft,pawprtvol,hdc,weiss)
78 :
79 : use defs_basis
80 : use m_errors
81 : use m_abicore
82 :
83 : use m_crystal, only : crystal_t
84 : use m_green, only : green_type,init_green,destroy_green
85 : use m_paw_dmft, only : paw_dmft_type
86 : use m_oper, only : oper_type,init_oper,destroy_oper,print_oper
87 : use m_matlu, only : matlu_type,sym_matlu, print_matlu, gather_matlu,&
88 : & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,copy_matlu,slm2ylm_matlu
89 : use m_hu, only : destroy_vee,hu_type,init_vee,rotatevee_hu,vee_type
90 : use m_datafordmft, only : compute_levels
91 :
92 : !Arguments ------------------------------------
93 : !scalars
94 : ! type(pawang_type), intent(in) :: pawang
95 : type(crystal_t),intent(in) :: cryst_struc
96 : type(green_type), intent(inout) :: green
97 : type(paw_dmft_type), intent(in) :: paw_dmft
98 : type(hu_type), intent(inout) :: hu(cryst_struc%ntypat)
99 : type(oper_type), intent(inout) :: hdc
100 : integer, intent(in) :: pawprtvol
101 : type(green_type), intent(inout) :: weiss
102 :
103 : !Local variables ------------------------------
104 : type :: level2_type
105 : integer, pointer :: repart(:,:) => null()
106 : integer, ABI_CONTIGUOUS pointer :: ocp(:,:) => null()
107 : integer, ABI_CONTIGUOUS pointer :: transition(:,:) => null()
108 : integer, ABI_CONTIGUOUS pointer :: transition_m(:,:) => null()
109 : end type level2_type
110 : type :: level1_type
111 : real(dp), pointer :: config(:) => null()
112 : end type level1_type
113 : ! scalars
114 : character(len=500) :: message
115 : integer :: iatom,ifreq,im,im1,isppol,ispinor,ispinor1
116 : integer :: lpawu,mbandc,natom,ndim,nkpt,nspinor,nsppol,nsppol_imp,testblock,useylm
117 : ! complex(dp) :: g,g0,w
118 : ! arrays
119 69 : complex(dp), allocatable :: Id(:,:,:,:)
120 69 : type(matlu_type), allocatable :: eigvectmatlu(:)
121 69 : type(matlu_type), allocatable :: udens_atoms(:)
122 69 : type(oper_type) :: energy_level
123 69 : type(green_type) :: green_hubbard
124 69 : type(matlu_type), allocatable :: level_diag(:)
125 : complex(dp) :: omega_current
126 69 : type(vee_type), allocatable :: vee_rotated(:)
127 : ! ************************************************************************
128 69 : mbandc=paw_dmft%mbandc
129 69 : nkpt=paw_dmft%nkpt
130 69 : nsppol=paw_dmft%nsppol
131 69 : natom=paw_dmft%natom
132 69 : nspinor=paw_dmft%nspinor
133 :
134 :
135 : !Initialise for compiler
136 69 : omega_current=czero
137 :
138 : !======================================
139 : !Allocations: levels and eigenvectors
140 : !======================================
141 509 : ABI_MALLOC(level_diag,(natom))
142 440 : ABI_MALLOC(eigvectmatlu,(natom))
143 440 : ABI_MALLOC(udens_atoms,(natom))
144 69 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag)
145 69 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,eigvectmatlu)
146 69 : call init_matlu(natom,2,1,paw_dmft%lpawu,udens_atoms)
147 371 : do iatom=1,cryst_struc%natom
148 302 : lpawu=paw_dmft%lpawu(iatom)
149 371 : if(lpawu/=-1) then
150 4301 : level_diag(iatom)%mat=czero
151 : end if
152 : end do
153 :
154 69 : call init_oper(paw_dmft,energy_level,opt_ksloc=3)
155 69 : call compute_levels(energy_level,hdc,paw_dmft)
156 : !!========================
157 : !!Get KS eigenvalues
158 : !!========================
159 : ! call init_oper(paw_dmft,energy_level,opt_ksloc=3)
160 : ! do iband=1,mbandc
161 : ! do ikpt=1,nkpt
162 : ! do isppol=1,nsppol
163 : !! Take \epsilon_{nks}
164 : !! ========================
165 : ! energy_level%ks(isppol,ikpt,iband,iband)=paw_dmft%eigen_dft(isppol,ikpt,iband)
166 : ! end do
167 : ! end do
168 : ! end do
169 : !
170 : !
171 : !!======================================================================
172 : !!Compute atomic levels from projection of \epsilon_{nks} and symetrize
173 : !!======================================================================
174 : ! call loc_oper(energy_level,paw_dmft,1)
175 : ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels before sym and only DFT"
176 : ! call wrtout(std_out,message,'COLL')
177 : ! call print_matlu(energy_level%matlu,natom,1)
178 : ! do iatom = 1 , natom
179 : ! lpawu=paw_dmft%lpawu(iatom)
180 : ! if(lpawu/=-1) then
181 : ! do isppol=1,nsppol
182 : ! do ispinor=1,nspinor
183 : ! do im1=1,2*lpawu+1
184 : ! energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)=&
185 : !& energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)&
186 : !& -hdc%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)-paw_dmft%fermie
187 : ! end do
188 : ! end do
189 : ! end do
190 : !! write(std_out,*) "DC,fermie",hdc%matlu(iatom)%mat(1,1,1,1,1),paw_dmft%fermie
191 : ! end if
192 : ! end do ! natom
193 : ! call sym_matlu(cryst_struc,energy_level%matlu,pawang)
194 : !
195 : ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels for Fermi Level=",paw_dmft%fermie
196 : ! call wrtout(std_out,message,'COLL')
197 : !!call print_oper(energy_level,1,paw_dmft,1)
198 : ! call print_matlu(energy_level%matlu,natom,1)
199 :
200 : !========================
201 : !Compute Weiss function
202 : !========================
203 276 : ABI_MALLOC(Id,(20,20,nspinor,nspinor))
204 371 : do iatom = 1 , natom
205 302 : lpawu=paw_dmft%lpawu(iatom)
206 371 : if(lpawu/=-1) then
207 41867 : Id=czero
208 486 : do im=1,2*lpawu+1
209 893 : do ispinor=1,nspinor
210 800 : Id(im,im,ispinor,ispinor)=cone
211 : end do
212 : end do ! ib
213 93 : ndim = 2*lpawu+1
214 5925 : do ifreq=1,weiss%nw
215 5832 : if(weiss%w_type=="imag") then
216 5832 : omega_current=cmplx(zero,weiss%omega(ifreq),kind=dp)
217 0 : else if(green%w_type=="real") then
218 0 : omega_current=cmplx(weiss%omega(ifreq),zero,kind=dp)
219 : end if
220 27613 : do im=1,2*lpawu+1
221 115656 : do im1=1,2*lpawu+1
222 240296 : do isppol=1,nsppol
223 358880 : do ispinor=1,nspinor
224 430616 : do ispinor1=1,nspinor
225 : weiss%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)=&
226 : & ( omega_current*Id(im,im1,ispinor,ispinor1) - &
227 300144 : & energy_level%matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol))
228 : end do ! ispinor1
229 : end do ! ispinor
230 : end do ! isppol
231 : end do ! im1
232 : end do ! im
233 : end do ! ifreq
234 : end if ! lpawu
235 : end do ! natom
236 69 : ABI_FREE(Id)
237 :
238 : !=================================================================
239 : !Diagonalizes atomic levels and keep eigenvectors in eigvectmatlu
240 : !=================================================================
241 : !if jpawu=0, rotatevee_hu will have no effect so it is not necessary to
242 : !have a single rotation matrix for up and dn spins.
243 :
244 69 : if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2
245 69 : if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1
246 : ! Diagonalize energy levels
247 69 : useylm=paw_dmft%dmft_blockdiag
248 69 : if(useylm==1) call slm2ylm_matlu(energy_level%matlu,natom,paw_dmft,1,pawprtvol)
249 69 : testblock=1
250 0 : if(useylm==1) testblock=8
251 : call diag_matlu(energy_level%matlu,level_diag,natom,prtopt=pawprtvol,eigvectmatlu=eigvectmatlu, &
252 69 : & nsppol_imp=nsppol_imp,opt_real=1,test=testblock)
253 :
254 : ! Use rotation matrix to rotate interaction
255 440 : ABI_MALLOC(vee_rotated,(natom))
256 69 : call init_vee(paw_dmft,vee_rotated)
257 69 : if(useylm==1) then
258 0 : call rotatevee_hu(hu,paw_dmft,pawprtvol,eigvectmatlu,4,udens_atoms,vee_rotated)
259 : else
260 69 : call rotatevee_hu(hu,paw_dmft,pawprtvol,eigvectmatlu,1,udens_atoms,vee_rotated)
261 : endif
262 69 : call destroy_vee(paw_dmft,vee_rotated)
263 371 : ABI_FREE(vee_rotated)
264 : !write(std_out,*)"udens after rotatevee", udens_atoms(1)%value
265 69 : write(message,'(a,2x,a,f13.5)') ch10,&
266 138 : & " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie
267 69 : call wrtout(std_out,message,'COLL')
268 69 : call print_matlu(level_diag,natom,1)
269 :
270 : ! Print out
271 69 : if(nspinor==2) then
272 2 : write(message,'(a,2x,a,f13.5)') ch10,&
273 4 : & " == Print weiss for small freq"
274 2 : call wrtout(std_out,message,'COLL')
275 2 : call print_matlu(weiss%oper(1)%matlu,natom,1)
276 2 : write(message,'(a,2x,a,f13.5)') ch10,&
277 4 : & " == Print weiss for large freq"
278 2 : call wrtout(std_out,message,'COLL')
279 2 : call print_matlu(weiss%oper(weiss%nw)%matlu,natom,1)
280 : end if
281 :
282 : !========================
283 : !Compute Green function
284 : !========================
285 69 : call init_green(green_hubbard,paw_dmft,opt_oper_ksloc=2,wtype=green%w_type) ! initialize only matlu
286 : !write(std_out,*)"udens", udens_atoms(1)%value
287 : ! write(std_out,*)"levels", level_diag(1)%mat
288 69 : call green_atomic_hubbard(cryst_struc,green_hubbard,hu,level_diag,paw_dmft,udens_atoms)
289 : !call rotate_matlu(energy_level%matlu,natom,pawprtvol=3)
290 : !write(81,*) "I1",paw_dmft%omega_lo(1), real(green%oper(1)%matlu(1)%mat(1,1,1,1,1)),imag(green%oper(1)%matlu(1)%mat(1,1,1,1,1))
291 : !========================================================================
292 : !Rotate back Green function in the original basis before diagonalization
293 : !========================================================================
294 : !call print_matlu(level_diag,natom,1)
295 : !test scall rotate_matlu(level_diag,eigvectmatlu,natom,3)
296 : !todo_ab: add check here for back rotation
297 : !call print_matlu(level_diag,natom,1)
298 69 : write(message,'(2a,f13.5)') ch10," == Green function before rotation"
299 69 : call wrtout(std_out,message,'COLL')
300 69 : call print_matlu(green_hubbard%oper(1)%matlu,natom,1)
301 5133 : do ifreq=1,green_hubbard%nw
302 5064 : call rotate_matlu(green_hubbard%oper(ifreq)%matlu,eigvectmatlu,natom,0)
303 5064 : if(useylm==1) call slm2ylm_matlu(green_hubbard%oper(ifreq)%matlu,natom,paw_dmft,2,0)
304 5133 : call copy_matlu(green_hubbard%oper(ifreq)%matlu,green%oper(ifreq)%matlu,natom)
305 : end do
306 69 : write(message,'(2a,f13.5)') ch10," == Green function after rotation"
307 69 : call wrtout(std_out,message,'COLL')
308 69 : call print_matlu(green%oper(1)%matlu,natom,1)
309 69 : if(nspinor==2) then
310 2 : write(message,'(a,2x,a,f13.5)') ch10,&
311 4 : & " == Print green for small freq"
312 2 : call wrtout(std_out,message,'COLL')
313 2 : call print_matlu(green%oper(1)%matlu,natom,1)
314 2 : write(message,'(a,2x,a,f13.5)') ch10,&
315 4 : & " == Print green for large freq"
316 2 : call wrtout(std_out,message,'COLL')
317 2 : call print_matlu(green%oper(green%nw)%matlu,natom,1)
318 : end if
319 : !do ifreq=1,paw_dmft%dmft_nwlo
320 : !g=green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
321 : !g0=cone/weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
322 : !w=cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)
323 : !write(160,*) paw_dmft%omega_lo(ifreq),real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),imag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
324 : !write(161,*) paw_dmft%omega_lo(ifreq),real(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) ),imag(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) )
325 : !write(164,*) paw_dmft%omega_lo(ifreq),real(one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) ),imag(one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
326 : !write(166,*) paw_dmft%omega_lo(ifreq),real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)-one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) ),imag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)-one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
327 : !write(167,*) paw_dmft%omega_lo(ifreq),real((g-g0)/(g0*g)),imag((g-g0)/(g0*g))
328 : !write(168,*) paw_dmft%omega_lo(ifreq),real(1/g-w),imag(1/g-w)
329 : !write(169,*) paw_dmft%omega_lo(ifreq),real(1/g0-w),imag(1/g0-w)
330 : !write(170,*) paw_dmft%omega_lo(ifreq),w
331 : !write(171,*) paw_dmft%omega_lo(ifreq),real(1/g),imag(1/g)
332 : !write(172,*) paw_dmft%omega_lo(ifreq),real(w),imag(w)
333 :
334 : !! voir si en faisant GG0/(G-G0) cela reduit l'erreur
335 : !enddo
336 : !call abi_abort('COLL')
337 :
338 :
339 : !write(message,'(2a,f13.5)') ch10," == Print Energy levels after diagonalisation"
340 : !call wrtout(std_out,message,'COLL')
341 : !call print_matlu(energy_level%matlu,natom,1)
342 :
343 : !======================================
344 : !Deallocations and destroys
345 : !======================================
346 69 : call destroy_green(green_hubbard)
347 69 : call destroy_oper(energy_level)
348 69 : call destroy_matlu(level_diag,natom)
349 69 : call destroy_matlu(udens_atoms,natom)
350 69 : call destroy_matlu(eigvectmatlu,natom)
351 371 : ABI_FREE(level_diag)
352 371 : ABI_FREE(eigvectmatlu)
353 371 : ABI_FREE(udens_atoms)
354 : !!***
355 :
356 : contains
357 :
358 : !!****f* m_hubbard_one/green_atomic_hubbard
359 : !! NAME
360 : !! green_atomic_hubbard
361 : !!
362 : !! FUNCTION
363 : !!
364 : !!
365 : !! COPYRIGHT
366 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
367 : !! This file is distributed under the terms of the
368 : !! GNU General Public License, see ~abinit/COPYING
369 : !! or http://www.gnu.org/copyleft/gpl.txt .
370 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
371 : !!
372 : !! INPUTS
373 : !! cryst_struc
374 : !! istep = step of iteration for DFT.
375 : !! dft_occup
376 : !! mpi_enreg=information about MPI parallelization
377 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
378 : !! pawang <type(pawang)>=paw angular mesh and related data
379 : !!
380 : !! OUTPUT
381 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
382 : !!
383 : !! NOTES
384 : !!
385 : !! SOURCE
386 :
387 138 : subroutine green_atomic_hubbard(cryst_struc,green_hubbard,hu,level_diag,paw_dmft,udens_atoms)
388 :
389 :
390 69 : use defs_basis
391 : use m_errors
392 : use m_abicore
393 : use m_crystal, only : crystal_t
394 : use m_special_funcs, only : factorial, permutations
395 : use m_green, only : green_type,init_green,destroy_green
396 : use m_hu, only : hu_type
397 : use m_paw_dmft, only : paw_dmft_type
398 :
399 : !Arguments ------------------------------------
400 : !scalars
401 : ! type(pawang_type), intent(in) :: pawang
402 : type(crystal_t),intent(in) :: cryst_struc
403 : type(green_type), intent(inout) :: green_hubbard
404 : type(paw_dmft_type), intent(in) :: paw_dmft
405 : type(matlu_type),intent(in) :: level_diag(cryst_struc%natom)
406 : type(hu_type), intent(inout) :: hu(cryst_struc%ntypat)
407 : type(matlu_type),intent(in) :: udens_atoms(cryst_struc%natom)
408 :
409 : !Local variables ------------------------------
410 : ! scalars
411 : integer :: cnk,iacc,iatom,iconfig,ielec,ifreq,ilevel,im,im1,isppol,ispinor,itrans,jconfig,jelec
412 : integer :: lpawu,m_temp,nconfig,ndim,nelec,nlevels,nspinor,nsppol,occupied_level,sum_test,gpu_option
413 : integer :: nconfig_next
414 69 : integer, allocatable :: occup(:,:),nconfig_nelec(:)
415 : character(len=500) :: message
416 69 : integer, ABI_CONTIGUOUS pointer :: ocp(:,:),ocp_next(:,:),transition(:,:),transition_m(:,:)
417 : ! arrays
418 69 : type(green_type) :: green_hubbard_realw
419 69 : type(level2_type), allocatable :: occ_level(:)
420 69 : type(level1_type), allocatable :: e_nelec(:)
421 69 : complex(dp), allocatable :: green_temp(:,:)
422 69 : complex(dp), allocatable :: green_temp_realw(:,:)
423 : complex(dp) :: Z_part
424 69 : real(dp), allocatable :: maxener(:),minener(:)
425 69 : real(dp), allocatable :: elevels(:)
426 : real(dp) :: emax,emin,eshift,prtopt, Ej_np1, Ei_n,beta,maxarg_exp,tmp
427 : !************************************************************************
428 : ABI_NVTX_START_RANGE(NVTX_DMFT_HUBBARD_ONE)
429 69 : maxarg_exp=300
430 :
431 : ! hu is not used anymore.
432 : if(hu(1)%lpawu==0) then
433 : end if
434 : ! ======================================
435 : ! General loop over atoms
436 : ! ======================================
437 69 : nsppol=paw_dmft%nsppol
438 69 : nspinor=paw_dmft%nspinor
439 69 : prtopt=1
440 69 : beta=one/paw_dmft%temp
441 69 : call init_green(green_hubbard_realw,paw_dmft,opt_oper_ksloc=2,wtype=green_hubbard%w_type) ! initialize only matlu
442 69 : gpu_option=paw_dmft%gpu_option
443 :
444 371 : do iatom=1,cryst_struc%natom
445 302 : lpawu=paw_dmft%lpawu(iatom)
446 371 : if(lpawu/=-1) then
447 93 : nlevels=nsppol*nspinor*(2*lpawu+1)
448 93 : if(nsppol==1.and.nspinor==1) nlevels=2*nspinor*(2*lpawu+1)
449 :
450 : ! ===================================
451 : ! Allocations
452 : ! ===================================
453 1158 : ABI_MALLOC(occ_level,(0:nlevels))
454 279 : ABI_MALLOC(maxener,(0:nlevels))
455 186 : ABI_MALLOC(minener,(0:nlevels))
456 279 : ABI_MALLOC(elevels,(nlevels))
457 1158 : ABI_MALLOC(e_nelec,(0:nlevels))
458 972 : do nelec=0,nlevels ! number of electrons
459 879 : cnk=nint(permutations(nlevels,nelec)/factorial(nelec))
460 3516 : ABI_MALLOC(occ_level(nelec)%repart ,(cnk,nelec))
461 3516 : ABI_MALLOC(occ_level(nelec)%ocp ,(cnk,nlevels))
462 3516 : ABI_MALLOC(occ_level(nelec)%transition ,(cnk,nlevels-nelec))
463 2637 : ABI_MALLOC(occ_level(nelec)%transition_m,(cnk,nlevels))
464 2637 : ABI_MALLOC(e_nelec (nelec)%config ,(cnk))
465 119244 : e_nelec(nelec)%config(:)=zero
466 : ! write(std_out,*) "permutations",nint(permutations(nlevels,nelec)/factorial(nelec))
467 : ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,1)
468 : ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,2)
469 : ! for a given nb of electrons nelec, gives for a given repartition
470 : ! of electron, the position of the ielec electron inside atomic
471 : ! levels
472 : ! levels
473 : end do
474 372 : ABI_MALLOC(occup,(0:nlevels,nlevels))
475 279 : ABI_MALLOC(nconfig_nelec,(0:nlevels))
476 :
477 : ! ===================================
478 : ! Initialization
479 : ! ===================================
480 972 : nconfig_nelec=0
481 93 : nconfig=1
482 8789 : occup=0
483 93 : nconfig_nelec(0)=1
484 879 : occup(0,:)=0
485 879 : iacc=0
486 879 : elevels=zero
487 : ndim=2*lpawu+1
488 237 : do isppol=1,nsppol
489 383 : do ispinor=1,nspinor
490 956 : do im1=1,(2*lpawu+1)
491 666 : iacc=iacc+1
492 666 : elevels(iacc)=level_diag(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor-1)*ndim,isppol)
493 666 : if(abs(aimag(level_diag(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor-1)*ndim,isppol)))>tol8) then
494 0 : message = " Hubbard I: levels are imaginary"
495 0 : write(std_out,*) level_diag(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor-1)*ndim,isppol)
496 0 : ABI_BUG(message)
497 : end if
498 812 : if(nsppol==1.and.nspinor==1) then
499 120 : elevels(nspinor*(2*lpawu+1)+iacc)=level_diag(iatom)%mat(im1+(ispinor-1)*ndim,im1+(ispinor-1)*ndim,isppol)
500 : end if
501 : ! ! change it by: real(level_diag) with warning
502 : end do
503 : end do
504 : end do
505 :
506 : ! ===================================
507 : ! Compute possible occupations
508 : ! ===================================
509 : ! Value for nelec=0:
510 : nconfig_nelec(0)=1
511 879 : occ_level(0)%ocp(1,:)=0
512 : ! Loop on possible occupation of levels with nelec
513 879 : do nelec=1,nlevels ! number of electrons
514 : ! write(message,'(2a,i3,a)') ch10," For number of electrons", &
515 : ! & nelec," positions of electrons are:"
516 : ! call wrtout(std_out,message,'COLL')
517 : ! write(std_out,*) "nelec",nelec
518 : ! write(std_out,*) "nlevels",nlevels
519 786 : call combin(1,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup)
520 786 : if(nconfig_nelec(nelec)/=nint(permutations(nlevels,nelec)/factorial(nelec))) then
521 0 : message = " BUG in hubbard_one/combin"
522 0 : ABI_BUG(message)
523 : end if
524 1441748 : occ_level(nelec)%ocp=zero
525 119058 : do iconfig=1,nconfig_nelec(nelec)
526 836277 : do ielec=1,nelec
527 : ! occ_level%repart: gives the place of electron ielec for the configuration iconfig (among the config for the total number of electron nelec
528 717312 : occupied_level=occ_level(nelec)%repart(iconfig,ielec)
529 : ! occ_level%ocp: gives if level occupied_level is occupied or not
530 835491 : occ_level(nelec)%ocp(iconfig,occupied_level)=1
531 : end do
532 : end do
533 : end do
534 :
535 : ! ===================================
536 : ! Print possible occupations
537 : ! ===================================
538 : if(prtopt>3) then
539 : do nelec=0,nlevels ! number of electrons f
540 : write(message,'(2a,i3,2a,i5,a)') ch10," For",nelec," electrons, ", &
541 : & "there are ",nconfig_nelec(nelec)," repartitions which are:"
542 : call wrtout(std_out,message,'COLL')
543 : do iconfig=1,nconfig_nelec(nelec)
544 : write(message,'(40i4)') (occ_level(nelec)%ocp(iconfig,ilevel),ilevel=1,nlevels),&
545 : & (occ_level(nelec)%repart(iconfig,ielec),ielec=1,nelec)
546 : call wrtout(std_out,message,'COLL')
547 : end do
548 : end do
549 : end if
550 :
551 : ! ============================================
552 : ! Compute energy for each of the occupations
553 : ! ============================================
554 972 : do nelec=0,nlevels !
555 119151 : e_nelec(nelec)%config=zero
556 119151 : do iconfig=1,nconfig_nelec(nelec)
557 : ! First compute energy level contribution
558 835584 : do ielec=1,nelec
559 : e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) &
560 835584 : & + elevels(occ_level(nelec)%repart(iconfig,ielec))
561 : end do
562 : ! write(std_out,*) "Nelec",nelec,"iconfig",iconfig,"eleve",e_nelec(nelec)%config(iconfig)
563 :
564 : ! Second: Compute interaction part
565 : ! do ielec=1,nelec-1 ! compute interaction among the nelec electrons in the configuration iconfig
566 : ! e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) &
567 : ! & + hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), &
568 : ! & occ_level(nelec)%repart(iconfig,ielec+1))
569 : ! enddo
570 836463 : do ielec=1,nelec ! compute interaction among the nelec electrons in the configuration iconfig
571 5682944 : do jelec=1,nelec
572 : e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) &
573 : ! & + hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), &
574 : & + dble(udens_atoms(iatom)%mat(occ_level(nelec)%repart(iconfig,ielec), &
575 5564672 : & occ_level(nelec)%repart(iconfig,jelec),1))/2.d0 ! udens(i,i)=0
576 : ! write(std_out,*) ielec,occ_level(nelec)%repart(iconfig,ielec)
577 : ! write(std_out,*) jelec,occ_level(nelec)%repart(iconfig,jelec)
578 : ! write(std_out,*)hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), &
579 : ! & occ_level(nelec)%repart(iconfig,jelec))/2.d0
580 : end do ! jelec
581 : end do ! ielec
582 : ! write(std_out,*) "Nelec",nelec,"iconfig",iconfig,"ecorr",e_nelec(nelec)%config(iconfig)
583 :
584 : end do ! iconfig
585 120030 : maxener(nelec)=maxval(-e_nelec(nelec)%config(:))
586 120123 : minener(nelec)=minval(-e_nelec(nelec)%config(:))
587 : end do
588 : ! write(std_out,*) "maxener", maxener(:)
589 1065 : emax=maxval(maxener(:))
590 : emin=minval(minener(:))
591 : eshift=zero
592 93 : eshift=emax/two
593 93 : eshift=emax-maxarg_exp/beta
594 : ! eshift=emax
595 : ! write(std_out,*)"emax",emax
596 : ! write(std_out,*)"emin",emin
597 : ! write(std_out,*)"eshift",eshift
598 93 : write(message,'(a,3x,3a,3x,a)') ch10," Hubbard I: Energies as a", &
599 93 : & " function of number of electrons",ch10,&
600 186 : & " Nelec Min. Ene. Max. Ener."
601 93 : call wrtout(std_out,message,'COLL')
602 972 : do nelec=0,nlevels
603 879 : write(message,'(3x,a,i4,2f17.7)') "HI", nelec,&
604 240060 : & minval(e_nelec(nelec)%config(:)),maxval(e_nelec(nelec)%config(:))
605 972 : call wrtout(std_out,message,'COLL')
606 : end do
607 :
608 : ! ===================================
609 : ! Print possibles occupations
610 : ! ===================================
611 : if(prtopt>3) then
612 : do nelec=0,nlevels ! number of electrons
613 : write(message,'(2a,i3,2a,i5,3a)') ch10," For",nelec," electrons, ", &
614 : & "there are ",nconfig_nelec(nelec)," repartitions which are :", &
615 : & ch10,"Energy and Occupations"
616 : call wrtout(std_out,message,'COLL')
617 : do iconfig=1,nconfig_nelec(nelec)
618 : write(message,'(f12.6,20i4)') e_nelec(nelec)%config(iconfig),&
619 : & (occ_level(nelec)%repart(iconfig,ielec),ielec=1,nelec)
620 : call wrtout(std_out,message,'COLL')
621 : end do
622 : end do
623 : end if
624 :
625 : ! sum_test=zero
626 : ! do ielec=1,nelec+1
627 : ! sum_test = sum_test + (occ_level(nelec)%repart(iconfig,ielec) &
628 : ! & -occ_level(nelec)%repart(iconfig,ielec))
629 : ! enddo
630 : ! ===================================
631 : ! Built transitions between configurations
632 : ! ===================================
633 879 : do nelec=0,nlevels-1
634 :
635 :
636 :
637 786 : nconfig = nconfig_nelec(nelec)
638 786 : nconfig_next = nconfig_nelec(nelec+1)
639 786 : transition => occ_level(nelec)%transition
640 786 : transition_m => occ_level(nelec)%transition_m
641 786 : ocp => occ_level(nelec)%ocp
642 786 : ocp_next => occ_level(nelec+1)%ocp
643 :
644 786 : if(gpu_option==ABI_GPU_DISABLED) then
645 : !$OMP PARALLEL DO PRIVATE(iconfig,jconfig,itrans,sum_test,m_temp)
646 118965 : do iconfig=1,nconfig
647 : itrans=0 ! transition from iconfig
648 158149325 : do jconfig=1, nconfig_next
649 158030360 : sum_test=0
650 2337281800 : do ilevel=1,nlevels
651 : ! test if their is one electron added to the starting configuration
652 : sum_test=sum_test + &
653 2179251440 : & (ocp_next(jconfig,ilevel)-ocp(iconfig,ilevel))**2
654 : ! save the level for the electron added
655 2337281800 : if(ocp_next(jconfig,ilevel)==1.and.ocp(iconfig,ilevel)==0) then
656 2179251440 : m_temp=ilevel
657 : end if
658 : end do ! ilevel
659 158148539 : if(sum_test==1) then
660 717312 : itrans=itrans+1
661 : !if(itrans>nlevels-nelec) then
662 : ! write(message,'(a,4i4)') "BUG: itrans is to big in hubbard_one",itrans,iconfig,jconfig,ilevel
663 : ! call wrtout(std_out,message,'COLL')
664 : !end if
665 717312 : transition(iconfig,itrans)=jconfig ! jconfig=config(n+1) obtained after transition
666 717312 : transition_m(iconfig,itrans)=m_temp ! level to fill to do the transition
667 : end if
668 : end do ! jconfig
669 : end do ! iconfig
670 : else if(gpu_option==ABI_GPU_OPENMP) then
671 : #ifdef HAVE_OPENMP_OFFLOAD
672 : !$OMP TARGET PARALLEL DO PRIVATE(iconfig,jconfig,itrans,sum_test,m_temp) &
673 : !$OMP MAP(tofrom:transition,transition_m) MAP(to:ocp,ocp_next)
674 : do iconfig=1,nconfig
675 : itrans=0 ! transition from iconfig
676 : do jconfig=1, nconfig_next
677 : sum_test=0
678 : do ilevel=1,nlevels
679 : ! test if their is one electron added to the starting configuration
680 : sum_test=sum_test + &
681 : & (ocp_next(jconfig,ilevel)-ocp(iconfig,ilevel))**2
682 : ! save the level for the electron added
683 : if(ocp_next(jconfig,ilevel)==1.and.ocp(iconfig,ilevel)==0) then
684 : m_temp=ilevel
685 : end if
686 : end do ! ilevel
687 : if(sum_test==1) then
688 : itrans=itrans+1
689 : !if(itrans>nlevels-nelec) then
690 : ! write(message,'(a,4i4)') "BUG: itrans is to big in hubbard_one",itrans,iconfig,jconfig,ilevel
691 : ! call wrtout(std_out,message,'COLL')
692 : !end if
693 : transition(iconfig,itrans)=jconfig ! jconfig=config(n+1) obtained after transition
694 : transition_m(iconfig,itrans)=m_temp ! level to fill to do the transition
695 : end if
696 : end do ! jconfig
697 : end do ! iconfig
698 : #endif
699 : end if
700 :
701 :
702 :
703 93 : if(prtopt>3) then
704 : do iconfig=1,nconfig
705 : write(std_out,'(a,2i5,a,18i5)') "occ_level", nelec,&
706 : & iconfig," :",(transition(iconfig,itrans),itrans=1,nlevels-nelec)
707 : write(std_out,'(a,2i5,a,18i5)') "electron added", nelec,iconfig,&
708 : & " :",(transition_m(iconfig,itrans),itrans=1,nlevels-nelec)
709 : end do
710 : end if
711 :
712 : end do ! nelec
713 :
714 : ! ===================================
715 : ! Built Partition Function
716 : ! ===================================
717 93 : Z_part=czero
718 : ! do nelec=1,nlevels-1
719 972 : do nelec=0,nlevels
720 119244 : do iconfig=1,nconfig_nelec(nelec)
721 118272 : Ei_n = e_nelec (nelec )%config(iconfig) + eshift
722 119151 : Z_part=Z_part+dexp(-Ei_n*beta)
723 : ! write(std_out,*) "fonction de partition",nelec,iconfig, Z_part,Ei_n*beta,Ei_n,eshift
724 : end do
725 : end do
726 : ! write(std_out,*) "Z_part",Z_part
727 :
728 : ! ===================================
729 : ! Built Green Function
730 : ! ===================================
731 372 : ABI_MALLOC(green_temp,(green_hubbard%nw,nlevels))
732 279 : ABI_MALLOC(green_temp_realw,(green_hubbard%nw,nlevels))
733 : ! For each freq.
734 :
735 44255 : green_temp=czero
736 44255 : green_temp_realw=czero
737 93 : tmp=zero
738 879 : do nelec=0,nlevels-1
739 : ! write(std_out,*) "For nelec =",nelec
740 119058 : do iconfig=1,nconfig_nelec(nelec)
741 : ! write(std_out,*) "The config nb:",iconfig
742 836277 : do itrans=1,nlevels-nelec
743 717312 : jconfig = occ_level(nelec )%transition(iconfig,itrans)
744 717312 : m_temp = occ_level(nelec )%transition_m(iconfig,itrans)
745 717312 : Ej_np1 = e_nelec (nelec+1)%config(jconfig) + eshift
746 717312 : Ei_n = e_nelec (nelec )%config(iconfig) + eshift
747 : ! write(std_out,'(a,i4,a)') "Transition nb:",itrans,"involve"
748 : ! write(std_out,'(a,i4,a)') " jconfig=",jconfig
749 : ! write(std_out,'(a,i4,a)') " m_temp=",m_temp
750 39909283 : do ifreq=1,green_hubbard%nw
751 39073792 : if(green_hubbard%w_type=="imag") then
752 39073792 : omega_current=cmplx(zero,green_hubbard%omega(ifreq),kind=dp)
753 0 : else if(green_hubbard%w_type=="real") then
754 0 : omega_current=cmplx(green_hubbard%omega(ifreq),zero,kind=dp)
755 : end if
756 : green_temp(ifreq,m_temp)=green_temp(ifreq,m_temp)+ &
757 : & (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))/ &
758 39073792 : & ( omega_current +Ei_n-Ej_np1)
759 : if(ifreq==1.and.m_temp==1) tmp=tmp+(dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))
760 :
761 : green_temp_realw(ifreq,m_temp)=green_temp_realw(ifreq,m_temp)+ &
762 : & (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))/ &
763 39791104 : & ( omega_current +Ei_n-Ej_np1)
764 : end do
765 : ! green_temp_realw(m_temp)=green_temp_realw(m_temp)+ &
766 : ! & (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta)) -> will give one at the end
767 : ! write(std_out,*) "green",-Ej_np1*beta,-Ei_n*beta,dexp(-Ej_np1*beta),dexp(-Ei_n*beta)
768 : end do
769 : end do
770 : end do
771 : ! write(std_out,*) "tmp",tmp
772 : ilevel=0
773 188 : do ispinor=1,nspinor
774 334 : do isppol=1,nsppol
775 907 : do im=1,(2*lpawu+1)
776 666 : ilevel=ilevel+1
777 : ! write(std_out,'(16e15.6)') paw_dmft%omega_lo(ifreq),(real(green_temp_realw(ilevel)/Z_part),ilevel=1,nlevels)
778 32188 : do ifreq=1,green_hubbard%nw
779 31376 : green_hubbard%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)=green_temp(ifreq,ilevel)/Z_part
780 32042 : green_hubbard_realw%oper(ifreq)%matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)=green_temp_realw(ifreq,ilevel)/Z_part
781 : end do
782 : end do
783 : end do
784 : end do
785 :
786 : ! End calculation for this frequency
787 93 : ABI_FREE(green_temp)
788 93 : ABI_FREE(green_temp_realw)
789 :
790 : ! ===================================
791 : ! Deallocations
792 : ! ===================================
793 972 : do nelec=0,nlevels
794 879 : ABI_FREE(occ_level(nelec)%repart)
795 879 : ABI_FREE(occ_level(nelec)%ocp)
796 879 : ABI_FREE(occ_level(nelec)%transition)
797 879 : ABI_FREE(occ_level(nelec)%transition_m)
798 972 : ABI_FREE(e_nelec(nelec)%config)
799 : end do
800 93 : ABI_FREE(occ_level)
801 93 : ABI_FREE(occup)
802 93 : ABI_FREE(nconfig_nelec)
803 93 : ABI_FREE(e_nelec)
804 93 : ABI_FREE(elevels)
805 93 : ABI_FREE(maxener)
806 93 : ABI_FREE(minener)
807 : end if
808 : end do
809 69 : call destroy_green(green_hubbard_realw)
810 :
811 : ABI_NVTX_END_RANGE()
812 :
813 138 : end subroutine green_atomic_hubbard
814 : !!***
815 :
816 : !!****f* m_hubbard_one/combin
817 : !! NAME
818 : !! combin
819 : !!
820 : !! FUNCTION
821 : !!
822 : !!
823 : !! COPYRIGHT
824 : !! Copyright (C) 1999-2026 ABINIT group (BAmadon)
825 : !! This file is distributed under the terms of the
826 : !! GNU General Public License, see ~abinit/COPYING
827 : !! or http://www.gnu.org/copyleft/gpl.txt .
828 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
829 : !!
830 : !! INPUTS
831 : !!
832 : !!
833 : !! OUTPUT
834 : !!
835 : !! NOTES
836 : !!
837 : !! SOURCE
838 :
839 118179 : recursive subroutine combin(ielec,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup)
840 :
841 : !Arguments ------------------------------------
842 : !scalars
843 : ! type(pawang_type), intent(in) :: pawang
844 : integer, intent(in) :: ielec,nelec,nlevels
845 : integer, intent(inout) :: nconfig,nconfig_nelec(0:nlevels)
846 : integer, intent(inout) :: occup(0:nlevels,nlevels)
847 : ! type :: level2_type
848 : ! integer, pointer :: repart(:,:)
849 : ! end type
850 : type(level2_type), intent(inout) :: occ_level(0:nlevels)
851 : ! integer, intent(in) :: prtopt
852 :
853 : !Local variables ------------------------------
854 : ! scalars
855 : integer :: max_ielec,pos,min_ielec,jelec,prtopt
856 : character(len=500) :: message
857 : ! arrays
858 : !************************************************************************
859 118179 : prtopt=1
860 118179 : max_ielec=nlevels-nelec+ielec
861 : ! write(std_out,*) "call to combin ielec,nelec,nlevels",ielec,nelec,nlevels
862 : select case (ielec)
863 : case (1)
864 117393 : min_ielec=1
865 : case default
866 118179 : min_ielec=occup(nelec,ielec-1)+1
867 : end select
868 : ! write(std_out,*) "For ielec", ielec, "min_ielec,max_ielec",min_ielec,max_ielec
869 353751 : do pos = min_ielec, max_ielec
870 353751 : if(ielec==nelec) then
871 118179 : occup(nelec,ielec)=pos
872 118179 : nconfig=nconfig+1
873 118179 : nconfig_nelec(nelec)=nconfig_nelec(nelec)+1
874 : ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,1)
875 : ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,2)
876 835491 : do jelec=1,nelec
877 : ! write(std_out,*) "nconfig",nconfig_nelec(nelec),nelec
878 : ! write(std_out,*) "occup",occup(nelec,jelec)
879 835491 : occ_level(nelec)%repart(nconfig_nelec(nelec),jelec)=occup(nelec,jelec)
880 : end do
881 : ! write(std_out,*) "For ielec", ielec, "case nelec"
882 : if(prtopt>=3) then
883 : write(message,'(a,i3,a,30i5)') "For ielec",ielec," Occupf are", (occup(nelec,jelec),jelec=1,nelec)
884 : call wrtout(std_out,message,'COLL')
885 : end if
886 : else
887 117393 : occup(nelec,ielec)=pos
888 : ! write(std_out,*) "For ielec", ielec, "case 1 and default"
889 117393 : call combin(ielec+1,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup)
890 : if(prtopt>=3) then
891 : write(message,'(a,i3,a,30i5)') "For ielec",ielec," Occup are", (occup(nelec,jelec),jelec=1,nelec)
892 : call wrtout(std_out,message,'COLL')
893 : end if
894 : end if
895 : end do
896 :
897 69 : end subroutine combin
898 : !!***
899 :
900 : end subroutine hubbard_one
901 : !!***
902 :
903 : END MODULE m_hubbard_one
904 : !!***
|
