TABLE OF CONTENTS
ABINIT/hubbard_one [ Functions ]
NAME
hubbard_one
FUNCTION
Solve the hubbard one approximation
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (BAmadon) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cryst_struc istep = step of iteration for LDA. lda_occup mpi_enreg=informations about MPI parallelization paw_dmft = data for self-consistent LDA+DMFT calculations. pawang <type(pawang)>=paw angular mesh and related data
OUTPUT
paw_dmft = data for self-consistent LDA+DMFT calculations.
NOTES
PARENTS
impurity_solve,spectral_function
CHILDREN
combin,destroy_green,init_green,wrtout
SOURCE
37 #if defined HAVE_CONFIG_H 38 #include "config.h" 39 #endif 40 41 42 #include "abi_common.h" 43 44 subroutine hubbard_one(cryst_struc,green,hu,paw_dmft,pawang,pawprtvol,hdc,weiss) 45 46 use defs_basis 47 use m_errors 48 use m_profiling_abi 49 50 use m_pawang, only : pawang_type 51 use m_crystal, only : crystal_t 52 use m_green, only : green_type,init_green,destroy_green 53 use m_paw_dmft, only : paw_dmft_type 54 use m_oper, only : oper_type,init_oper,destroy_oper,loc_oper,print_oper 55 use m_matlu, only : matlu_type,sym_matlu, print_matlu, gather_matlu,& 56 & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,copy_matlu 57 use m_hu, only : hu_type,rotatevee_hu 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'hubbard_one' 63 use interfaces_14_hidewrite 64 use interfaces_68_dmft, except_this_one => hubbard_one 65 !End of the abilint section 66 67 implicit none 68 69 !Arguments ------------------------------------ 70 !scalars 71 ! type(pawang_type), intent(in) :: pawang 72 type(crystal_t),intent(in) :: cryst_struc 73 type(green_type), intent(inout) :: green 74 type(paw_dmft_type), intent(in) :: paw_dmft 75 type(hu_type), intent(inout) :: hu(cryst_struc%ntypat) 76 type(pawang_type), intent(in) :: pawang 77 type(oper_type), intent(inout) :: hdc 78 integer, intent(in) :: pawprtvol 79 type(green_type), intent(inout) :: weiss 80 81 !Local variables ------------------------------ 82 type :: level2_type 83 integer, pointer :: repart(:,:) => null() 84 integer, pointer :: ocp(:,:) => null() 85 integer, pointer :: transition(:,:) => null() 86 integer, pointer :: transition_m(:,:) => null() 87 end type level2_type 88 type :: level1_type 89 real(dp), pointer :: config(:) => null() 90 end type level1_type 91 ! scalars 92 character(len=500) :: message 93 integer :: iatom,ifreq,im,im1,isppol,ispinor,ispinor1 94 integer :: lpawu,mbandc,natom,nkpt,nspinor,nsppol,nsppol_imp,tndim 95 ! complex(dpc) :: g,g0,w 96 ! arrays 97 complex(dp), allocatable :: Id(:,:,:,:) 98 type(coeff2c_type), allocatable :: eigvectmatlu(:,:) 99 type(coeff2_type), allocatable :: udens_atoms(:) 100 type(oper_type) :: energy_level 101 type(green_type) :: green_hubbard 102 type(matlu_type), allocatable :: level_diag(:) 103 complex(dpc) :: omega_current 104 ! ************************************************************************ 105 mbandc=paw_dmft%mbandc 106 nkpt=paw_dmft%nkpt 107 nsppol=paw_dmft%nsppol 108 natom=paw_dmft%natom 109 nspinor=paw_dmft%nspinor 110 111 112 !Initialise for compiler 113 omega_current=czero 114 115 !====================================== 116 !Allocations: levels and eigenvectors 117 !====================================== 118 ABI_DATATYPE_ALLOCATE(level_diag,(natom)) 119 ABI_DATATYPE_ALLOCATE(eigvectmatlu,(natom,nsppol)) 120 ABI_DATATYPE_ALLOCATE(udens_atoms,(natom)) 121 call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag) 122 do iatom=1,cryst_struc%natom 123 lpawu=paw_dmft%lpawu(iatom) 124 if(lpawu/=-1) then 125 tndim=nspinor*(2*lpawu+1) 126 do isppol=1,nsppol 127 ABI_ALLOCATE(eigvectmatlu(iatom,isppol)%value,(tndim,tndim)) 128 end do 129 ABI_ALLOCATE(udens_atoms(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1))) 130 level_diag(iatom)%mat=czero 131 end if 132 end do 133 134 call init_oper(paw_dmft,energy_level,opt_ksloc=3) 135 call compute_levels(cryst_struc,energy_level,hdc,pawang,paw_dmft) 136 !!======================== 137 !!Get KS eigenvalues 138 !!======================== 139 ! call init_oper(paw_dmft,energy_level,opt_ksloc=3) 140 ! do iband=1,mbandc 141 ! do ikpt=1,nkpt 142 ! do isppol=1,nsppol 143 !! Take \epsilon_{nks} 144 !! ======================== 145 ! energy_level%ks(isppol,ikpt,iband,iband)=paw_dmft%eigen_lda(isppol,ikpt,iband) 146 ! end do 147 ! end do 148 ! end do 149 ! 150 ! 151 !!====================================================================== 152 !!Compute atomic levels from projection of \epsilon_{nks} and symetrize 153 !!====================================================================== 154 ! call loc_oper(energy_level,paw_dmft,1) 155 ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels before sym and only LDA" 156 ! call wrtout(std_out,message,'COLL') 157 ! call print_matlu(energy_level%matlu,natom,1) 158 ! do iatom = 1 , natom 159 ! lpawu=paw_dmft%lpawu(iatom) 160 ! if(lpawu/=-1) then 161 ! do isppol=1,nsppol 162 ! do ispinor=1,nspinor 163 ! do im1=1,2*lpawu+1 164 ! energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)=& 165 !& energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)& 166 !& -hdc%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)-paw_dmft%fermie 167 ! end do 168 ! end do 169 ! end do 170 !! write(std_out,*) "DC,fermie",hdc%matlu(iatom)%mat(1,1,1,1,1),paw_dmft%fermie 171 ! end if 172 ! end do ! natom 173 ! call sym_matlu(cryst_struc,energy_level%matlu,pawang) 174 ! 175 ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels for Fermi Level=",paw_dmft%fermie 176 ! call wrtout(std_out,message,'COLL') 177 !!call print_oper(energy_level,1,paw_dmft,1) 178 ! call print_matlu(energy_level%matlu,natom,1) 179 180 !======================== 181 !Compute Weiss function 182 !======================== 183 ABI_ALLOCATE(Id,(20,20,nspinor,nspinor)) 184 do iatom = 1 , natom 185 lpawu=paw_dmft%lpawu(iatom) 186 if(lpawu/=-1) then 187 Id=czero 188 do im=1,2*lpawu+1 189 do ispinor=1,nspinor 190 Id(im,im,ispinor,ispinor)=cone 191 end do 192 end do ! ib 193 do ifreq=1,weiss%nw 194 if(weiss%w_type=="imag") then 195 omega_current=cmplx(zero,weiss%omega(ifreq),kind=dp) 196 else if(green%w_type=="real") then 197 omega_current=cmplx(weiss%omega(ifreq),zero,kind=dp) 198 end if 199 do im=1,2*lpawu+1 200 do im1=1,2*lpawu+1 201 do isppol=1,nsppol 202 do ispinor=1,nspinor 203 do ispinor1=1,nspinor 204 weiss%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=& 205 & ( omega_current*Id(im,im1,ispinor,ispinor1) - & 206 & energy_level%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)) 207 end do ! ispinor1 208 end do ! ispinor 209 end do ! isppol 210 end do ! im1 211 end do ! im 212 end do ! ifreq 213 end if ! lpawu 214 end do ! natom 215 ABI_DEALLOCATE(Id) 216 217 !================================================================= 218 !Diagonalizes atomic levels and keep eigenvectors in eigvectmatlu 219 !================================================================= 220 !if jpawu=0, rotatevee_hu will have no effect so it is not necessary to 221 !have a single rotation matrix for up and dn spins. 222 223 if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2 224 if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1 225 ! Diagonalize energy levels 226 call diag_matlu(energy_level%matlu,level_diag,natom,& 227 & prtopt=pawprtvol,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,optreal=1) 228 229 ! Use rotation matrix to rotate interaction 230 call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,1) 231 !write(std_out,*)"udens after rotatevee", udens_atoms(1)%value 232 write(message,'(a,2x,a,f13.5)') ch10,& 233 & " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie 234 call wrtout(std_out,message,'COLL') 235 call print_matlu(level_diag,natom,1) 236 237 ! Print out 238 if(nspinor==2) then 239 write(message,'(a,2x,a,f13.5)') ch10,& 240 & " == Print weiss for small freq" 241 call wrtout(std_out,message,'COLL') 242 call print_matlu(weiss%oper(1)%matlu,natom,1) 243 write(message,'(a,2x,a,f13.5)') ch10,& 244 & " == Print weiss for large freq" 245 call wrtout(std_out,message,'COLL') 246 call print_matlu(weiss%oper(weiss%nw)%matlu,natom,1) 247 end if 248 249 !======================== 250 !Compute Green function 251 !======================== 252 call init_green(green_hubbard,paw_dmft,opt_oper_ksloc=2,wtype=green%w_type) ! initialize only matlu 253 !write(std_out,*)"udens", udens_atoms(1)%value 254 ! write(std_out,*)"levels", level_diag(1)%mat 255 call green_atomic_hubbard(cryst_struc,green_hubbard,hu,level_diag,paw_dmft,udens_atoms) 256 !call rotate_matlu(energy_level%matlu,natom,pawprtvol=3) 257 !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)) 258 !======================================================================== 259 !Rotate back Green function in the original basis before diagonalization 260 !======================================================================== 261 !call print_matlu(level_diag,natom,1) 262 !test scall rotate_matlu(level_diag,eigvectmatlu,natom,3) 263 !todo_ab: add check here for back rotation 264 !call print_matlu(level_diag,natom,1) 265 write(message,'(2a,f13.5)') ch10," == Green function before rotation" 266 call wrtout(std_out,message,'COLL') 267 call print_matlu(green_hubbard%oper(1)%matlu,natom,1) 268 do ifreq=1,green_hubbard%nw 269 call rotate_matlu(green_hubbard%oper(ifreq)%matlu,eigvectmatlu,natom,3,0) 270 call copy_matlu(green_hubbard%oper(ifreq)%matlu,green%oper(ifreq)%matlu,natom) 271 end do 272 write(message,'(2a,f13.5)') ch10," == Green function after rotation" 273 call wrtout(std_out,message,'COLL') 274 call print_matlu(green%oper(1)%matlu,natom,1) 275 if(nspinor==2) then 276 write(message,'(a,2x,a,f13.5)') ch10,& 277 & " == Print green for small freq" 278 call wrtout(std_out,message,'COLL') 279 call print_matlu(green%oper(1)%matlu,natom,1) 280 write(message,'(a,2x,a,f13.5)') ch10,& 281 & " == Print green for large freq" 282 call wrtout(std_out,message,'COLL') 283 call print_matlu(green%oper(green%nw)%matlu,natom,1) 284 end if 285 !do ifreq=1,paw_dmft%dmft_nwlo 286 !g=green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) 287 !g0=cone/weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) 288 !w=cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp) 289 !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)) 290 !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) ) 291 !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)) 292 !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)) 293 !write(167,*) paw_dmft%omega_lo(ifreq),real((g-g0)/(g0*g)),imag((g-g0)/(g0*g)) 294 !write(168,*) paw_dmft%omega_lo(ifreq),real(1/g-w),imag(1/g-w) 295 !write(169,*) paw_dmft%omega_lo(ifreq),real(1/g0-w),imag(1/g0-w) 296 !write(170,*) paw_dmft%omega_lo(ifreq),w 297 !write(171,*) paw_dmft%omega_lo(ifreq),real(1/g),imag(1/g) 298 !write(172,*) paw_dmft%omega_lo(ifreq),real(w),imag(w) 299 300 !! voir si en faisant GG0/(G-G0) cela reduit l'erreur 301 !enddo 302 !call leave_new('COLL') 303 304 305 !write(message,'(2a,f13.5)') ch10," == Print Energy levels after diagonalisation" 306 !call wrtout(std_out,message,'COLL') 307 !call print_matlu(energy_level%matlu,natom,1) 308 309 !====================================== 310 !Deallocations and destroys 311 !====================================== 312 call destroy_green(green_hubbard) 313 call destroy_oper(energy_level) 314 call destroy_matlu(level_diag,natom) 315 ABI_DATATYPE_DEALLOCATE(level_diag) 316 do iatom=1,cryst_struc%natom 317 lpawu=paw_dmft%lpawu(iatom) 318 if(lpawu/=-1) then 319 ABI_DEALLOCATE(udens_atoms(iatom)%value) 320 do isppol=1,nsppol 321 ABI_DEALLOCATE(eigvectmatlu(iatom,isppol)%value) 322 end do 323 end if 324 end do 325 ABI_DATATYPE_DEALLOCATE(eigvectmatlu) 326 ABI_DATATYPE_DEALLOCATE(udens_atoms) 327 328 CONTAINS
hubbard_one/combin [ Functions ]
[ Top ] [ hubbard_one ] [ Functions ]
NAME
combin
FUNCTION
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (BAmadon) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
OUTPUT
NOTES
PARENTS
CHILDREN
SOURCE
777 recursive subroutine combin(ielec,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup) 778 779 use defs_basis 780 781 !This section has been created automatically by the script Abilint (TD). 782 !Do not modify the following lines by hand. 783 #undef ABI_FUNC 784 #define ABI_FUNC 'combin' 785 use interfaces_14_hidewrite 786 !End of the abilint section 787 788 implicit none 789 790 !Arguments ------------------------------------ 791 !scalars 792 ! type(pawang_type), intent(in) :: pawang 793 integer, intent(in) :: ielec,nelec,nlevels 794 integer, intent(inout) :: nconfig,nconfig_nelec(0:nlevels) 795 integer, intent(inout) :: occup(0:nlevels,nlevels) 796 ! type :: level2_type 797 ! integer, pointer :: repart(:,:) 798 ! end type 799 type(level2_type), intent(inout) :: occ_level(0:nlevels) 800 ! integer, intent(in) :: prtopt 801 802 !Local variables ------------------------------ 803 ! scalars 804 integer :: max_ielec,pos,min_ielec,jelec,prtopt 805 character(len=500) :: message 806 ! arrays 807 !************************************************************************ 808 prtopt=1 809 max_ielec=nlevels-nelec+ielec 810 ! write(std_out,*) "call to combin ielec,nelec,nlevels",ielec,nelec,nlevels 811 select case (ielec) 812 case (1) 813 min_ielec=1 814 case default 815 min_ielec=occup(nelec,ielec-1)+1 816 end select 817 ! write(std_out,*) "For ielec", ielec, "min_ielec,max_ielec",min_ielec,max_ielec 818 do pos = min_ielec, max_ielec 819 if(ielec==nelec) then 820 occup(nelec,ielec)=pos 821 nconfig=nconfig+1 822 nconfig_nelec(nelec)=nconfig_nelec(nelec)+1 823 ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,1) 824 ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,2) 825 do jelec=1,nelec 826 ! write(std_out,*) "nconfig",nconfig_nelec(nelec),nelec 827 ! write(std_out,*) "occup",occup(nelec,jelec) 828 occ_level(nelec)%repart(nconfig_nelec(nelec),jelec)=occup(nelec,jelec) 829 end do 830 ! write(std_out,*) "For ielec", ielec, "case nelec" 831 if(prtopt>=3) then 832 write(message,'(a,i3,a,30i5)') "For ielec",ielec," Occupf are", (occup(nelec,jelec),jelec=1,nelec) 833 call wrtout(std_out,message,'COLL') 834 end if 835 else 836 occup(nelec,ielec)=pos 837 ! write(std_out,*) "For ielec", ielec, "case 1 and default" 838 call combin(ielec+1,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup) 839 if(prtopt>=3) then 840 write(message,'(a,i3,a,30i5)') "For ielec",ielec," Occup are", (occup(nelec,jelec),jelec=1,nelec) 841 call wrtout(std_out,message,'COLL') 842 end if 843 end if 844 end do 845 846 end subroutine combin 847 848 end subroutine hubbard_one
hubbard_one/green_atomic_hubbard [ Functions ]
[ Top ] [ hubbard_one ] [ Functions ]
NAME
green_atomic_hubbard
FUNCTION
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (BAmadon) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cryst_struc istep = step of iteration for LDA. lda_occup mpi_enreg=informations about MPI parallelization paw_dmft = data for self-consistent LDA+DMFT calculations. pawang <type(pawang)>=paw angular mesh and related data
OUTPUT
paw_dmft = data for self-consistent LDA+DMFT calculations.
NOTES
PARENTS
hubbard_one
CHILDREN
combin,destroy_green,init_green,wrtout
SOURCE
366 subroutine green_atomic_hubbard(cryst_struc,green_hubbard,hu,level_diag,paw_dmft,udens_atoms) 367 368 369 use defs_basis 370 use m_errors 371 use m_profiling_abi 372 use m_crystal, only : crystal_t 373 use m_special_funcs, only : factorial, permutations 374 use m_green, only : green_type,init_green,destroy_green 375 use m_hu, only : hu_type 376 use m_paw_dmft, only : paw_dmft_type 377 378 !This section has been created automatically by the script Abilint (TD). 379 !Do not modify the following lines by hand. 380 #undef ABI_FUNC 381 #define ABI_FUNC 'green_atomic_hubbard' 382 use interfaces_14_hidewrite 383 !End of the abilint section 384 385 implicit none 386 387 !Arguments ------------------------------------ 388 !scalars 389 ! type(pawang_type), intent(in) :: pawang 390 type(crystal_t),intent(in) :: cryst_struc 391 type(green_type), intent(inout) :: green_hubbard 392 type(paw_dmft_type), intent(in) :: paw_dmft 393 type(matlu_type),intent(in) :: level_diag(cryst_struc%natom) 394 type(hu_type), intent(inout) :: hu(cryst_struc%ntypat) 395 type(coeff2_type),intent(in) :: udens_atoms(cryst_struc%natom) 396 397 !Local variables ------------------------------ 398 ! scalars 399 integer :: cnk,iacc,iatom,iconfig,ielec,ifreq,ilevel,im,im1,isppol,ispinor,itrans,jconfig,jelec 400 integer :: lpawu,m_temp,nconfig,nelec,nlevels,nspinor,nsppol,occupied_level,sum_test 401 integer, allocatable :: occup(:,:),nconfig_nelec(:) 402 character(len=500) :: message 403 ! arrays 404 type(green_type) :: green_hubbard_realw 405 type(level2_type), allocatable :: occ_level(:) 406 type(level1_type), allocatable :: e_nelec(:) 407 complex(dpc), allocatable :: green_temp(:,:) 408 complex(dpc), allocatable :: green_temp_realw(:,:) 409 complex(dpc) :: Z_part 410 real(dp), allocatable :: maxener(:),minener(:) 411 real(dp), allocatable :: elevels(:) 412 real(dp) :: emax,emin,eshift,prtopt, Ej_np1, Ei_n,beta,maxarg_exp,tmp 413 !************************************************************************ 414 maxarg_exp=300 415 416 ! hu is not used anymore. 417 if(hu(1)%lpawu==0) then 418 end if 419 ! ====================================== 420 ! General loop over atoms 421 ! ====================================== 422 nsppol=paw_dmft%nsppol 423 nspinor=paw_dmft%nspinor 424 prtopt=1 425 beta=one/paw_dmft%temp 426 call init_green(green_hubbard_realw,paw_dmft,opt_oper_ksloc=2,wtype=green_hubbard%w_type) ! initialize only matlu 427 428 do iatom=1,cryst_struc%natom 429 lpawu=paw_dmft%lpawu(iatom) 430 if(lpawu/=-1) then 431 nlevels=nsppol*nspinor*(2*lpawu+1) 432 if(nsppol==1.and.nspinor==1) nlevels=2*nspinor*(2*lpawu+1) 433 434 ! =================================== 435 ! Allocations 436 ! =================================== 437 ABI_DATATYPE_ALLOCATE(occ_level,(0:nlevels)) 438 ABI_ALLOCATE(maxener,(0:nlevels)) 439 ABI_ALLOCATE(minener,(0:nlevels)) 440 ABI_ALLOCATE(elevels,(nlevels)) 441 ABI_DATATYPE_ALLOCATE(e_nelec,(0:nlevels)) 442 do nelec=0,nlevels ! number of electrons 443 cnk=nint(permutations(nlevels,nelec)/factorial(nelec)) 444 ABI_ALLOCATE(occ_level(nelec)%repart ,(cnk,nelec)) 445 ABI_ALLOCATE(occ_level(nelec)%ocp ,(cnk,nlevels)) 446 ABI_ALLOCATE(occ_level(nelec)%transition ,(cnk,nlevels-nelec)) 447 ABI_ALLOCATE(occ_level(nelec)%transition_m,(cnk,nlevels)) 448 ABI_ALLOCATE(e_nelec (nelec)%config ,(cnk)) 449 e_nelec(nelec)%config(:)=zero 450 ! write(std_out,*) "permutations",nint(permutations(nlevels,nelec)/factorial(nelec)) 451 ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,1) 452 ! write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,2) 453 ! for a given nb of electrons nelec, gives for a given repartition 454 ! of electron, the position of the ielec electron inside atomic 455 ! levels 456 ! levels 457 end do 458 ABI_ALLOCATE(occup,(0:nlevels,nlevels)) 459 ABI_ALLOCATE(nconfig_nelec,(0:nlevels)) 460 461 ! =================================== 462 ! Initialization 463 ! =================================== 464 nconfig_nelec=0 465 nconfig=1 466 occup=0 467 nconfig_nelec(0)=1 468 occup(0,:)=0 469 iacc=0 470 elevels=zero 471 do isppol=1,nsppol 472 do ispinor=1,nspinor 473 do im1=1,(2*lpawu+1) 474 iacc=iacc+1 475 elevels(iacc)=level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor) 476 if(abs(aimag(level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor)))>tol8) then 477 message = " Hubbard I: levels are imaginary" 478 write(std_out,*) level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor) 479 MSG_BUG(message) 480 end if 481 if(nsppol==1.and.nspinor==1) then 482 elevels(nspinor*(2*lpawu+1)+iacc)=level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor) 483 end if 484 ! ! change it by: real(level_diag) with warning 485 end do 486 end do 487 end do 488 489 ! =================================== 490 ! Compute possible occupations 491 ! =================================== 492 ! Value for nelec=0: 493 nconfig_nelec(0)=1 494 occ_level(0)%ocp(1,:)=0 495 ! Loop on possible occupation of levels with nelec 496 do nelec=1,nlevels ! number of electrons 497 ! write(message,'(2a,i3,a)') ch10," For number of electrons", & 498 ! & nelec," positions of electrons are:" 499 ! call wrtout(std_out,message,'COLL') 500 ! write(std_out,*) "nelec",nelec 501 ! write(std_out,*) "nlevels",nlevels 502 call combin(1,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup) 503 if(nconfig_nelec(nelec)/=nint(permutations(nlevels,nelec)/factorial(nelec))) then 504 message = " BUG in hubbard_one/combin" 505 MSG_BUG(message) 506 end if 507 occ_level(nelec)%ocp=zero 508 do iconfig=1,nconfig_nelec(nelec) 509 do ielec=1,nelec 510 ! occ_level%repart: gives the place of electron ielec for the configuration iconfig (among the config for the total number of electron nelec 511 occupied_level=occ_level(nelec)%repart(iconfig,ielec) 512 ! occ_level%ocp: gives if level occupied_level is occupied or not 513 occ_level(nelec)%ocp(iconfig,occupied_level)=1 514 end do 515 end do 516 end do 517 518 ! =================================== 519 ! Print possible occupations 520 ! =================================== 521 if(prtopt>3) then 522 do nelec=0,nlevels ! number of electrons f 523 write(message,'(2a,i3,2a,i5,a)') ch10," For",nelec," electrons, ", & 524 & "there are ",nconfig_nelec(nelec)," repartitions which are:" 525 call wrtout(std_out,message,'COLL') 526 do iconfig=1,nconfig_nelec(nelec) 527 write(message,'(40i4)') (occ_level(nelec)%ocp(iconfig,ilevel),ilevel=1,nlevels),& 528 & (occ_level(nelec)%repart(iconfig,ielec),ielec=1,nelec) 529 call wrtout(std_out,message,'COLL') 530 end do 531 end do 532 end if 533 534 ! ============================================ 535 ! Compute energy for each of the occupations 536 ! ============================================ 537 do nelec=0,nlevels ! 538 e_nelec(nelec)%config=zero 539 do iconfig=1,nconfig_nelec(nelec) 540 ! First compute energy level contribution 541 do ielec=1,nelec 542 e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) & 543 & + elevels(occ_level(nelec)%repart(iconfig,ielec)) 544 end do 545 ! write(std_out,*) "Nelec",nelec,"iconfig",iconfig,"eleve",e_nelec(nelec)%config(iconfig) 546 547 ! Second: Compute interaction part 548 ! do ielec=1,nelec-1 ! compute interaction among the nelec electrons in the configuration iconfig 549 ! e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) & 550 ! & + hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), & 551 ! & occ_level(nelec)%repart(iconfig,ielec+1)) 552 ! enddo 553 do ielec=1,nelec ! compute interaction among the nelec electrons in the configuration iconfig 554 do jelec=1,nelec 555 e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) & 556 ! & + hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), & 557 & + udens_atoms(iatom)%value(occ_level(nelec)%repart(iconfig,ielec), & 558 & occ_level(nelec)%repart(iconfig,jelec))/2.d0 ! udens(i,i)=0 559 ! write(std_out,*) ielec,occ_level(nelec)%repart(iconfig,ielec) 560 ! write(std_out,*) jelec,occ_level(nelec)%repart(iconfig,jelec) 561 ! write(std_out,*)hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), & 562 ! & occ_level(nelec)%repart(iconfig,jelec))/2.d0 563 end do ! jelec 564 end do ! ielec 565 ! write(std_out,*) "Nelec",nelec,"iconfig",iconfig,"ecorr",e_nelec(nelec)%config(iconfig) 566 567 end do ! iconfig 568 maxener(nelec)=maxval(-e_nelec(nelec)%config(:)) 569 minener(nelec)=minval(-e_nelec(nelec)%config(:)) 570 end do 571 ! write(std_out,*) "maxener", maxener(:) 572 emax=maxval(maxener(:)) 573 emin=minval(minener(:)) 574 eshift=zero 575 eshift=emax/two 576 eshift=emax-maxarg_exp/beta 577 ! eshift=emax 578 ! write(std_out,*)"emax",emax 579 ! write(std_out,*)"emin",emin 580 ! write(std_out,*)"eshift",eshift 581 write(message,'(a,3x,3a,3x,a)') ch10," Hubbard I: Energies as a", & 582 & " function of number of electrons",ch10,& 583 & " Nelec Min. Ene. Max. Ener." 584 call wrtout(std_out,message,'COLL') 585 do nelec=0,nlevels 586 write(message,'(3x,a,i4,2f17.7)') "HI", nelec,& 587 & minval(e_nelec(nelec)%config(:)),maxval(e_nelec(nelec)%config(:)) 588 call wrtout(std_out,message,'COLL') 589 end do 590 591 ! =================================== 592 ! Print possibles occupations 593 ! =================================== 594 if(prtopt>3) then 595 do nelec=0,nlevels ! number of electrons 596 write(message,'(2a,i3,2a,i5,3a)') ch10," For",nelec," electrons, ", & 597 & "there are ",nconfig_nelec(nelec)," repartitions which are :", & 598 & ch10,"Energy and Occupations" 599 call wrtout(std_out,message,'COLL') 600 do iconfig=1,nconfig_nelec(nelec) 601 write(message,'(f12.6,20i4)') e_nelec(nelec)%config(iconfig),& 602 & (occ_level(nelec)%repart(iconfig,ielec),ielec=1,nelec) 603 call wrtout(std_out,message,'COLL') 604 end do 605 end do 606 end if 607 608 ! sum_test=zero 609 ! do ielec=1,nelec+1 610 ! sum_test = sum_test + (occ_level(nelec)%repart(iconfig,ielec) & 611 ! & -occ_level(nelec)%repart(iconfig,ielec)) 612 ! enddo 613 ! =================================== 614 ! Built transitions between configurations 615 ! =================================== 616 do nelec=0,nlevels-1 617 do iconfig=1,nconfig_nelec(nelec) 618 itrans=0 ! transition from iconfig 619 do jconfig=1, nconfig_nelec(nelec+1) 620 sum_test=0 621 do ilevel=1,nlevels 622 ! test if their is one electron added to the starting configuration 623 sum_test=sum_test + & 624 & (occ_level(nelec+1)%ocp(jconfig,ilevel)- occ_level(nelec)%ocp(iconfig,ilevel))**2 625 ! save the level for the electron added 626 if(occ_level(nelec+1)%ocp(jconfig,ilevel)==1.and.occ_level(nelec)%ocp(iconfig,ilevel)==0) then 627 m_temp=ilevel 628 end if 629 end do ! ilevel 630 if(sum_test==1) then 631 itrans=itrans+1 632 if(itrans>nlevels-nelec) then 633 write(message,'(a,4i4)') "BUG: itrans is to big in hubbard_one",itrans,iconfig,jconfig,ilevel 634 call wrtout(std_out,message,'COLL') 635 end if 636 occ_level(nelec)%transition(iconfig,itrans)=jconfig ! jconfig=config(n+1) obtained after transition 637 occ_level(nelec)%transition_m(iconfig,itrans)=m_temp ! level to fill to do the transition 638 end if 639 end do ! jconfig 640 if(prtopt>3) then 641 write(std_out,'(a,2i5,a,18i5)') "occ_level", nelec,& 642 & iconfig," :",(occ_level(nelec)%transition(iconfig,itrans),itrans=1,nlevels-nelec) 643 write(std_out,'(a,2i5,a,18i5)') "electron added", nelec,iconfig,& 644 & " :",(occ_level(nelec)%transition_m(iconfig,itrans),itrans=1,nlevels-nelec) 645 end if 646 end do ! iconfig 647 end do ! nelec 648 649 ! =================================== 650 ! Built Partition Function 651 ! =================================== 652 Z_part=czero 653 ! do nelec=1,nlevels-1 654 do nelec=0,nlevels 655 do iconfig=1,nconfig_nelec(nelec) 656 Ei_n = e_nelec (nelec )%config(iconfig) + eshift 657 Z_part=Z_part+dexp(-Ei_n*beta) 658 ! write(std_out,*) "fonction de partition",nelec,iconfig, Z_part,Ei_n*beta,Ei_n,eshift 659 end do 660 end do 661 ! write(std_out,*) "Z_part",Z_part 662 663 ! =================================== 664 ! Built Green Function 665 ! =================================== 666 ABI_ALLOCATE(green_temp,(green_hubbard%nw,nlevels)) 667 ABI_ALLOCATE(green_temp_realw,(green_hubbard%nw,nlevels)) 668 ! For each freq. 669 670 green_temp=czero 671 green_temp_realw=czero 672 tmp=zero 673 do nelec=0,nlevels-1 674 ! write(std_out,*) "For nelec =",nelec 675 do iconfig=1,nconfig_nelec(nelec) 676 ! write(std_out,*) "The config nb:",iconfig 677 do itrans=1,nlevels-nelec 678 jconfig = occ_level(nelec )%transition(iconfig,itrans) 679 m_temp = occ_level(nelec )%transition_m(iconfig,itrans) 680 Ej_np1 = e_nelec (nelec+1)%config(jconfig) + eshift 681 Ei_n = e_nelec (nelec )%config(iconfig) + eshift 682 ! write(std_out,'(a,i4,a)') "Transition nb:",itrans,"involve" 683 ! write(std_out,'(a,i4,a)') " jconfig=",jconfig 684 ! write(std_out,'(a,i4,a)') " m_temp=",m_temp 685 do ifreq=1,green_hubbard%nw 686 if(green_hubbard%w_type=="imag") then 687 omega_current=cmplx(zero,green_hubbard%omega(ifreq),kind=dp) 688 else if(green_hubbard%w_type=="real") then 689 omega_current=cmplx(green_hubbard%omega(ifreq),zero,kind=dp) 690 end if 691 green_temp(ifreq,m_temp)=green_temp(ifreq,m_temp)+ & 692 & (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))/ & 693 & ( omega_current +Ei_n-Ej_np1) 694 if(ifreq==1.and.m_temp==1) tmp=tmp+(dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta)) 695 696 green_temp_realw(ifreq,m_temp)=green_temp_realw(ifreq,m_temp)+ & 697 & (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))/ & 698 & ( omega_current +Ei_n-Ej_np1) 699 end do 700 ! green_temp_realw(m_temp)=green_temp_realw(m_temp)+ & 701 ! & (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta)) -> will give one at the end 702 ! write(std_out,*) "green",-Ej_np1*beta,-Ei_n*beta,dexp(-Ej_np1*beta),dexp(-Ei_n*beta) 703 end do 704 end do 705 end do 706 ! write(std_out,*) "tmp",tmp 707 ilevel=0 708 do ispinor=1,nspinor 709 do isppol=1,nsppol 710 do im=1,(2*lpawu+1) 711 ilevel=ilevel+1 712 ! write(std_out,'(16e15.6)') paw_dmft%omega_lo(ifreq),(real(green_temp_realw(ilevel)/Z_part),ilevel=1,nlevels) 713 do ifreq=1,green_hubbard%nw 714 green_hubbard%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=green_temp(ifreq,ilevel)/Z_part 715 green_hubbard_realw%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=green_temp_realw(ifreq,ilevel)/Z_part 716 end do 717 end do 718 end do 719 end do 720 721 ! End calculation for this frequency 722 ABI_DEALLOCATE(green_temp) 723 ABI_DEALLOCATE(green_temp_realw) 724 725 ! =================================== 726 ! Deallocations 727 ! =================================== 728 do nelec=0,nlevels 729 ABI_DEALLOCATE(occ_level(nelec)%repart) 730 ABI_DEALLOCATE(occ_level(nelec)%ocp) 731 ABI_DEALLOCATE(occ_level(nelec)%transition) 732 ABI_DEALLOCATE(occ_level(nelec)%transition_m) 733 ABI_DEALLOCATE(e_nelec(nelec)%config) 734 end do 735 ABI_DATATYPE_DEALLOCATE(occ_level) 736 ABI_DEALLOCATE(occup) 737 ABI_DEALLOCATE(nconfig_nelec) 738 ABI_DATATYPE_DEALLOCATE(e_nelec) 739 ABI_DEALLOCATE(elevels) 740 ABI_DEALLOCATE(maxener) 741 ABI_DEALLOCATE(minener) 742 end if 743 end do 744 call destroy_green(green_hubbard_realw) 745 746 end subroutine green_atomic_hubbard