Branch data 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, contiguous, pointer :: ocp(:,:) => null()
107 : : integer, contiguous, pointer :: transition(:,:) => null()
108 : : integer, 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, 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 : : !!***
|
