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