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