TABLE OF CONTENTS
- m_paw_correlations/calc_ubare
- m_paw_correlations/calc_vee
- m_paw_correlations/loc_orbmom_cal
- m_paw_correlations/m_paw_correlations
- m_paw_correlations/pawpuxinit
- m_paw_correlations/pawuenergy
- m_paw_correlations/pawxenergy
- m_paw_correlations/setnoccmmp
- m_paw_correlations/setrhoijpbe0
m_paw_correlations/calc_ubare [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
calc_ubare
FUNCTION
Calculate the bare interaction on atomic orbitals
INPUTS
itypatcor = value of itypat for correlated species lpawu = angular momentum for correlated species pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data: pawang %lmax=Maximum value of angular momentum l+1 %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data: %mesh_size=Dimension of radial mesh %rad(mesh_size)=The coordinates of all the points of the radial mesh %radfact(mesh_size)=Factor used to compute radial integrals
OUTPUT
NOTES
SOURCE
2767 subroutine calc_ubare(itypatcor,lpawu,pawang,pawrad,pawtab,rmax) 2768 2769 !Arguments ------------------------------------ 2770 integer, intent(in) :: itypatcor,lpawu 2771 type(pawang_type),intent(in) :: pawang 2772 type(pawrad_type),intent(in) :: pawrad 2773 type(pawtab_type),target,intent(in) :: pawtab 2774 real(dp), optional, intent(in) :: rmax 2775 2776 !Local variables ------------------------------ 2777 !scalars 2778 integer :: ilmn,ilmn1,iln,iln1,isel,isel1,itypat,jlmn,jlmn1,jln,jln1 2779 integer :: klm,klm1,klmn,klmn1,ll,lm0 2780 integer :: lmin,lmax,lmn2_size,mesh_size,meshsz,mm 2781 real(dp) :: norm,r_for_intg,rg,rg1,ubare,uint,uint_tmp 2782 character(len=800) :: message 2783 !arrays 2784 real(dp),allocatable :: ff(:),gg(:),phiphj(:),phiphj1(:) 2785 2786 !************************************************************************ 2787 2788 itypat=itypatcor 2789 write(message,'(11a,f12.4,2a,i7,2a,f12.4,2a,i7,2a,f12.4)') & 2790 & ch10," =======================================================================",ch10, & 2791 & " == Calculation of diagonal bare Coulomb interaction on ATOMIC orbitals ",ch10, & 2792 & " (it is assumed that the wavefunction for the first reference ",ch10, & 2793 & " energy in PAW atomic data is an atomic eigenvalue)",ch10,ch10, & 2794 & " Max value of the radius in atomic data file =", pawrad%rmax ,ch10, & 2795 & " Max value of the mesh in atomic data file =", pawrad%mesh_size,ch10, & 2796 & " PAW radius is =", pawtab%rpaw,ch10, & 2797 & " PAW value of the mesh for integration is =", pawrad%int_meshsz,ch10, & 2798 & " Integral of atomic wavefunction until rpaw =", pawtab%ph0phiint(1) 2799 if(.not.present(rmax)) then 2800 call wrtout(ab_out,message,'COLL') 2801 call wrtout(std_out,message,'COLL') 2802 end if 2803 2804 mesh_size=pawrad%mesh_size 2805 2806 ! Definition of the mesh used for integration. 2807 if(present(rmax)) then 2808 if(rmax>pawrad%rmax) then 2809 write(message, '(a)' ) 'calc_ubare: the radius cannot be larger than the maximum radius of the mesh' 2810 ABI_ERROR(message) 2811 end if 2812 meshsz=pawrad_ifromr(pawrad,rmax)+5 2813 r_for_intg=rmax 2814 else 2815 meshsz=pawtab%partialwave_mesh_size 2816 r_for_intg=pawrad%rad(meshsz) ! (we could use r_for_intg=-1) 2817 end if 2818 2819 lmn2_size=pawtab%lmn2_size 2820 ABI_MALLOC(ff,(mesh_size)) 2821 ABI_MALLOC(gg,(mesh_size)) 2822 ABI_MALLOC(phiphj,(mesh_size)) 2823 ABI_MALLOC(phiphj1,(mesh_size)) 2824 do klmn=1,lmn2_size 2825 ilmn=pawtab%indklmn(7,klmn);jlmn=pawtab%indklmn(8,klmn) 2826 ! Select lpawu and first projectors il=jl=lpawu and first proj only 2827 if (( pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. & 2828 & (-pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. & 2829 & (pawtab%indlmn(3,ilmn)==1).and.(pawtab%indlmn(3,jlmn)==1) ) then 2830 klm=pawtab%indklmn(1,klmn);iln=pawtab%indlmn(5,ilmn);jln=pawtab%indlmn(5,jlmn) 2831 lmin=pawtab%indklmn(3,klmn);lmax=pawtab%indklmn(4,klmn) 2832 phiphj(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln) 2833 !write(6,*) "A",klmn,pawtab%klmntomn(1,klmn),pawtab%klmntomn(2,klmn),& 2834 !&pawtab%indklmn(7,klmn),pawtab%indklmn(8,klmn),pawtab%klmntomn(3,klmn),pawtab%klmntomn(4,klmn) 2835 do ll=lmin,lmin,2 2836 lm0=ll*ll+ll+1 2837 ff(1:meshsz)=phiphj(1:meshsz) 2838 call simp_gen(norm,ff,pawrad,r_for_intg=r_for_intg) 2839 call poisson(ff,ll,pawrad,gg) 2840 do klmn1=klmn,lmn2_size 2841 ilmn1=pawtab%indklmn(7,klmn);jlmn1=pawtab%indklmn(8,klmn) 2842 ! Select lpawu and first projectors il=jl=lpawu and first proj only 2843 if (( pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. & 2844 & (-pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. & 2845 & (pawtab%indlmn(3,ilmn1)==1).and.(pawtab%indlmn(3,jlmn1)==1) ) then 2846 !write(6,*) "A1",klmn1,pawtab%klmntomn(1,klmn1),pawtab%klmntomn(2,klmn1),& 2847 !&pawtab%indklmn(7,klmn1),pawtab%indklmn(8,klmn1),pawtab%klmntomn(3,klmn1),pawtab%klmntomn(4,klmn1) 2848 klm1=pawtab%indklmn(1,klmn1);iln1=pawtab%indlmn(5,ilmn1);jln1=pawtab%indlmn(5,jlmn1) 2849 phiphj1(1:meshsz)=pawtab%phi(1:meshsz,iln1)*pawtab%phi(1:meshsz,jln1) 2850 uint_tmp=zero 2851 if ((ll==lmin)) then 2852 ff(1)=zero 2853 ff(2:meshsz)=phiphj1(2:meshsz)*gg(2:meshsz)*two/pawrad%rad(2:meshsz) 2854 call simp_gen(uint_tmp,ff,pawrad,r_for_intg=r_for_intg) 2855 end if 2856 uint=zero 2857 do mm=-ll,ll 2858 isel =pawang%gntselect(lm0+mm,klm) 2859 isel1=pawang%gntselect(lm0+mm,klm1) 2860 if (isel>0.and.isel1>0) then 2861 rg =pawang%realgnt(isel) 2862 rg1=pawang%realgnt(isel1) 2863 uint=uint+uint_tmp*rg*rg1*two_pi 2864 end if 2865 end do 2866 if((pawtab%indklmn(5,klmn)==pawtab%indklmn(6,klmn)).and.& 2867 & (pawtab%indklmn(5,klmn1)==pawtab%indklmn(6,klmn1)).and.& 2868 & (pawtab%indklmn(5,klmn)==pawtab%indklmn(5,klmn1))) then 2869 ubare=uint*Ha_eV 2870 end if 2871 end if 2872 end do 2873 end do 2874 end if 2875 end do 2876 ABI_FREE(gg) 2877 ABI_FREE(ff) 2878 ABI_FREE(phiphj) 2879 ABI_FREE(phiphj1) 2880 2881 write(message,'(a,3(a,f12.4,a),2a,f12.4,a)') ch10," For an atomic wfn truncated at rmax =",r_for_intg,ch10,& 2882 & " The norm of the wfn is =",norm,ch10,& 2883 & " The bare interaction (no renormalization) =",ubare," eV",ch10,& 2884 & " The bare interaction (for a renorm. wfn ) =",ubare/norm/norm," eV" 2885 call wrtout(ab_out,message,'COLL') 2886 call wrtout(std_out,message,'COLL') 2887 if(r_for_intg < 10_dp .and. .not.present(rmax)) then 2888 write(message,'(a,f6.2,4a)') ' ( WARNING: The radial mesh in the atomic data file is cut at',r_for_intg,ch10,& 2889 & ' Use XML atomic data files to compute the bare Coulomb interaction',ch10,& 2890 & ' on a true normalized atomic wavefunction )' 2891 call wrtout(ab_out,message,'COLL') 2892 call wrtout(std_out,message,'COLL') 2893 end if 2894 if(present(rmax)) then 2895 write(message,'(2a)') " =======================================================================",ch10 2896 call wrtout(ab_out,message,'COLL') 2897 call wrtout(std_out,message,'COLL') 2898 end if 2899 2900 end subroutine calc_ubare
m_paw_correlations/calc_vee [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
calc_vee
FUNCTION
Compute matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
(angular part computed from Gaunt coefficients)
INPUTS
jpawu= value of J lpawu= value of l on which DFT+U applies upawu= value of U
OUTPUT
vee(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
SOURCE
835 subroutine calc_vee(f4of2_sla,f6of2_sla,jpawu,lpawu,pawang,upawu,vee,prtvol) 836 837 !Arguments --------------------------------------------- 838 !scalars 839 integer,intent(in) :: lpawu 840 integer,optional,intent(in) :: prtvol 841 real(dp),intent(in) :: upawu,jpawu 842 real(dp),intent(inout) :: f4of2_sla,f6of2_sla 843 type(pawang_type), intent(in) :: pawang 844 !arrays 845 real(dp),intent(out) :: vee(2*lpawu+1,2*lpawu+1,2*lpawu+1,2*lpawu+1) 846 847 !Local variables --------------------------------------- 848 !scalars 849 integer :: isela,iselb 850 integer :: klm0u,klma,klmb,kyc,lkyc 851 integer :: lmkyc,lmpawu 852 integer :: m1,m11,m2,m21,m3,m31,m4,m41,prtvol_ 853 integer :: mkyc,sz1 854 real(dp) :: ak,f4of2,f6of2 855 character(len=500) :: message 856 !arrays 857 real(dp),allocatable :: fk(:) 858 859 ! ************************************************************************* 860 861 DBG_ENTER("COLL") 862 863 864 prtvol_ = 3 865 if (present(prtvol)) then 866 prtvol_ = prtvol 867 end if 868 ! Select only atoms with +U 869 if(lpawu/=-1) then 870 871 ! ====================================================================== 872 ! C-PAW+U: Matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]] 873 ! 1. angular part computed from Gaunt coefficients 874 ! -------------------------------------------------------------------- 875 ! a. compute F(k) 876 ! --------------------------------------------- 877 ABI_MALLOC(fk,(lpawu+1)) 878 fk(1)=upawu 879 ! cf Slater Physical Review 165, p 665 (1968) [[cite:Slater1958]] 880 ! write(std_out,*) "f4of2_sla",pawtab(itypat)%f4of2_sla 881 if(lpawu==0) then 882 fk(1)=fk(1) 883 else if(lpawu==1) then 884 fk(2)=jpawu*5._dp 885 else if(lpawu==2) then 886 ! f4of2=0._dp 887 if(f4of2_sla<-0.1_dp) then 888 f4of2=0.625_dp 889 f4of2_sla=f4of2 890 else 891 f4of2=f4of2_sla 892 end if 893 fk(2)=jpawu*14._dp/(One+f4of2) 894 fk(3)=fk(2)*f4of2 895 if(abs(prtvol_)>=2) then 896 write(message,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,& 897 & "Slater parameters F^0, F^2, F^4 are",fk(1),fk(2),fk(3) 898 call wrtout(std_out,message,'COLL') 899 end if 900 else if(lpawu==3) then 901 f4of2=0.6681_dp 902 f6of2=0.4943_dp 903 if(f4of2_sla<-0.1_dp) then 904 f4of2=0.6681_dp 905 f4of2_sla=f4of2 906 else 907 f4of2=f4of2_sla 908 end if 909 if(f6of2_sla<-0.1_dp) then 910 f6of2=0.4943_dp 911 f6of2_sla=f6of2 912 else 913 f6of2=f6of2_sla 914 end if 915 fk(2)=jpawu*6435._dp/(286._dp+195._dp*f4of2+250._dp*f6of2) 916 fk(3)=fk(2)*f4of2 917 fk(4)=fk(2)*f6of2 918 if(abs(prtvol_)>=2) then 919 write(std_out,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,& 920 & "Slater parameters F^0, F^2, F^4, F^6 are",fk(1),fk(2),fk(3),fk(4) 921 end if 922 else 923 write(message, '(a,i0,2a)' )& 924 & ' lpawu=',lpawu,ch10,& 925 & ' lpawu not equal to 0 ,1 ,2 or 3 is not allowed' 926 ABI_ERROR(message) 927 end if 928 929 ! b. Compute ak and vee. 930 ! --------------------------------------------- 931 ! if (allocated(vee)) then 932 ! ABI_DEALLOCATE(vee) 933 ! end if 934 sz1=2*lpawu+1 935 ! ABI_ALLOCATE(vee,(sz1,sz1,sz1,sz1)) 936 vee=zero 937 lmpawu=(lpawu-1)**2+2*(lpawu-1)+1 ! number of m value below correlated orbitals 938 klm0u=lmpawu*(lmpawu+1)/2 ! value of klmn just below correlated orbitals 939 ! --------- 4 loops for interaction matrix 940 do m1=-lpawu,lpawu 941 m11=m1+lpawu+1 942 do m2=-lpawu,m1 943 m21=m2+lpawu+1 944 ! klma= number of pair before correlated orbitals + 945 ! number of pair for m1 lower than correlated orbitals 946 ! (m1+lpawu+1)*(lpawu-1) + number of pairs for correlated orbitals 947 ! before (m1,m2) + number of pair for m2 lower than current value 948 klma=klm0u+m11*lmpawu+(m11-1)*m11/2+m21 949 do m3=-lpawu,lpawu 950 m31=m3+lpawu+1 951 do m4=-lpawu,m3 952 m41=m4+lpawu+1 953 klmb=klm0u+m31*lmpawu+(m31-1)*m31/2+m41 954 ! --------- loop on k=1,2,3 (4 if f orbitals) 955 do kyc=1,2*lpawu+1,2 956 lkyc=kyc-1 957 lmkyc=(lkyc+1)*(lkyc)+1 958 ak=zero 959 do mkyc=-lkyc,lkyc,1 960 isela=pawang%gntselect(lmkyc+mkyc,klma) 961 iselb=pawang%gntselect(lmkyc+mkyc,klmb) 962 if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb) 963 end do 964 ! ----- end loop on k=1,2,3 (4 if f orbitals) 965 ak=ak/(two*dble(lkyc)+one) 966 vee(m11,m31,m21,m41)=ak*fk(lkyc/2+1)+vee(m11,m31,m21,m41) 967 end do !kyc 968 vee(m11,m31,m21,m41)=vee(m11,m31,m21,m41)*four_pi 969 vee(m21,m31,m11,m41)=vee(m11,m31,m21,m41) 970 vee(m11,m41,m21,m31)=vee(m11,m31,m21,m41) 971 vee(m21,m41,m11,m31)=vee(m11,m31,m21,m41) 972 end do 973 end do 974 end do 975 end do 976 ABI_FREE(fk) 977 endif 978 979 end subroutine calc_vee
m_paw_correlations/loc_orbmom_cal [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
loc_orbmom_cal
FUNCTION
Calculate the orbital magnetic moments in PAW spheres
INPUTS
INPUTS
compute_dmat= flag: if 1, nocc_{m,mp} is computed
dimdmat=first dimension of dmatpawu array
dmatpawu(dimdmat,dimdmat,nsppol*nspinor,natpawu)=input density matrix to be copied into noccmpp
dmatudiag= flag controlling the use of diagonalization:
0: no diagonalization of nocc_{m,mp}
1: diagonalized nocc_{m,mp} matrix is printed
2: dmatpawu matrix is expressed in the basis where nocc_(m,mp} is diagonal
impose_dmat= flag: if 1, nocc_{m,mp} is replaced by dmatpawu
indsym(4,nsym,natom)=indirect indexing array for atom labels
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
comm_atom=--optional-- MPI communicator over atoms
my_natom=number of atoms treated by current processor
natom=number of atoms in cell
natpawu=number of atoms on which PAW+U is applied
nspinor=number of spinorial components of the wavefunctions
nsppol=number of independant spin components
nsym=number of symmetry elements in space group
ntypat=number of atom types
paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
pawang <type(pawang_type)>=paw angular mesh and related data
pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
spinat(3,matom)=initial spin of each atom, in unit of hbar/2
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
typat(natom)=type for each atom
useexexch=1 if local-exact-exchange is activated
usepawu= /=0 if PAW+U is activated
OUTPUT
printing the values of orbital magnetic moments for atoms in the output file
NOTES
PARENTS
m_outscfcv
CHILDREN
m_paw_correlations
SOURCE
2952 subroutine loc_orbmom_cal(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,& 2953 & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawrad,pawprtvol,pawrhoij,pawtab,& 2954 & spinat,symafm,typat,useexexch,usepawu, & 2955 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 2956 2957 !Arguments --------------------------------------------- 2958 !scalars 2959 integer,intent(in) :: compute_dmat,dimdmat,dmatudiag,impose_dmat,my_natom,natom,natpawu 2960 integer,intent(in) :: nspinor,nsppol,nsym,ntypat,useexexch,usepawu 2961 integer,optional,intent(in) :: comm_atom 2962 type(pawang_type),intent(in) :: pawang 2963 integer,intent(in) :: pawprtvol 2964 !arrays 2965 integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),typat(natom) 2966 integer,optional,target,intent(in) :: mpi_atmtab(:) 2967 real(dp),intent(in) :: dmatpawu(dimdmat,dimdmat,nspinor*nsppol,natpawu*impose_dmat) 2968 real(dp),intent(in) :: spinat(3,natom) 2969 type(paw_ij_type),intent(in) :: paw_ij(my_natom) 2970 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 2971 type(pawtab_type),intent(in) :: pawtab(ntypat) 2972 integer,pointer :: my_atmtab(:) 2973 !Local variables --------------------------------------- 2974 !scalars 2975 logical :: paral_atom,my_atmtab_allocated 2976 character(len=5) :: orb_char 2977 integer :: cplex_dij,im1,im2,ndij,itypat,my_comm_atom 2978 integer :: my_lcur,my_iatom,coor,isp,lmin,lmax,me_atom 2979 real(dp),allocatable :: my_l_occmat(:,:,:,:) 2980 complex(dpc),allocatable :: op_l(:,:,:),cmfoccmat(:,:,:) 2981 real(dp) :: orb_mom(3) 2982 real(dp) :: sum_orb_mom(3) 2983 complex(dpc) :: my_sls_val 2984 character(len=500) :: message 2985 type(paw_ij_type), ABI_CONTIGUOUS pointer :: paw_ij_all(:) 2986 type(pawrhoij_type),ABI_CONTIGUOUS pointer :: pawrhoij_all(:)
m_paw_correlations/m_paw_correlations [ Modules ]
NAME
m_paw_correlations
FUNCTION
This module contains several routines related to the treatment of electronic
correlations in the PAW approach (DFT+U, exact-exchange, ...).
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (BA,FJ,MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_paw_correlations 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_dtset 30 use m_linalg_interfaces 31 use m_special_funcs 32 33 use m_io_tools, only : open_file 34 use m_pawang, only : pawang_type,pawang_init,pawang_free 35 use m_pawrad, only : pawrad_type,simp_gen,nderiv_gen,pawrad_ifromr,poisson 36 use m_pawtab, only : pawtab_type,pawtab_nullify,pawtab_free,pawtab_set_flags 37 use m_pawrhoij, only : pawrhoij_type,pawrhoij_gather, pawrhoij_nullify, pawrhoij_free 38 use m_paw_ij, only : paw_ij_type,paw_ij_gather, paw_ij_free, paw_ij_nullify 39 use m_paw_sphharm, only : mat_mlms2jmj,mat_slm2ylm,slxyzs 40 use m_paw_io, only : pawio_print_ij 41 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 42 use m_copy, only : alloc_copy 43 44 implicit none 45 46 private 47 48 !public procedures. 49 public :: pawpuxinit ! Initialize some data for PAW+U/PAW+LocalExactExchange/PAW+DMFT 50 public :: calc_vee ! Compute vee for DFT+U 51 public :: pawuenergy ! Compute contributions to energy for PAW+U 52 public :: pawxenergy ! Compute contributions to energy for PAW+[local exact exchange] 53 public :: setnoccmmp ! Compute DFT+U density matrix nocc_{m,m_prime} or impose it 54 public :: setrhoijpbe0 ! Impose value of rhoij for using an auxiliairy file (PBE0 only) 55 public :: calc_ubare ! Calculate the bare interaction on atomic orbitals 56 public :: loc_orbmom_cal ! calculate local orbital magnetic moments 57 CONTAINS !========================================================================================
m_paw_correlations/pawpuxinit [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
pawpuxinit
FUNCTION
Initialize some starting values of several arrays used in PAW+U/+DMFT or local exact-exchange calculations A-define useful indices for DFT+U/local exact-exchange B-Compute overlap between atomic wavefunction C-Compute matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]] (angular part computed from Gaunt coefficients)
INPUTS
dmatpuopt= select expression for the density matrix exchmix= mixing factor for local exact-exchange is_dfpt=true if we are running a DFPT calculation jpawu(ntypat)= value of J llexexch(ntypat)= value of l on which local exact-exchange applies llpawu(ntypat)= value of l on which DFT+U applies ntypat=number of types of atoms in unit cell. pawang <type(pawang_type)>=paw angular mesh and related data %lmax=Maximum value of angular momentum l+1 %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients pawprtvol=output printing level for PAW pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data: upawu(ntypat)= value of U use_dmft = 0 no PAW+DMFT, =1 PAW+DMFT useexexch= 0 if no local exact-exchange; 1 if local exact-exchange usepawu= 0 if no DFT+U; /=0 if DFT+U
OUTPUT
pawtab <type(pawtab_type)>=paw tabulated data read at start: %euijkl=(3,lmn2_size,lmn2_size)= array for computing DFT+U terms without occupancies %ij_proj= nproj*(nproju+1)/2 %klmntomn(4,lmn2_size)= Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn) %lnproju(nproj)= value of ln for projectors on which paw+u/local exact-exchange acts. %nproju=number of projectors for orbitals on which paw+u/local exact-exchange acts. %phiphjint(pawtabitypat%ij_proj)=Integral of Phi(:,i)*Phi(:,j) for correlated orbitals. %usepawu=0 if no DFT+U; /=0 if DFT+U %useexexch=0 if no local exact-exchange; 1 if local exact-exchange === if usepawu/=0 %jpawu= value of J %upawu= value of U %vee(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals === if useexexch/=0 %fk %vex(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
SOURCE
115 subroutine pawpuxinit(dmatpuopt,exchmix,f4of2_sla,f6of2_sla,is_dfpt,jpawu,llexexch,llpawu,& 116 & nspinor,ntypat,option_interaction,pawang,pawprtvol,pawrad,pawtab,upawu,use_dmft,& 117 & useexexch,usepawu,& 118 & ucrpa,lmagCalc) ! optional argument 119 120 !Arguments --------------------------------------------- 121 !scalars 122 integer,intent(in) :: dmatpuopt,nspinor,ntypat,pawprtvol,use_dmft,useexexch,usepawu 123 !Option for interaction energy in case of non-collinear magnetism: 124 ! 1: E_int=-J/4.N.(N-2) 125 ! 2: E_int=-J/2.(Nup.(Nup-1)+Ndn.(Ndn-1)) (Nup and Ndn are ill-defined) 126 ! 3: E_int=-J/4.( N.(N-2) + mx^2 + my^2 + mz^2 ) 127 ! Default is 3 128 integer,intent(in) :: option_interaction 129 logical :: is_dfpt 130 real(dp),intent(in) :: exchmix 131 type(pawang_type), intent(in) :: pawang 132 integer,optional, intent(in) :: ucrpa 133 !arrays 134 integer,intent(in) :: llexexch(ntypat),llpawu(ntypat) 135 real(dp),intent(in) :: jpawu(ntypat),upawu(ntypat) 136 real(dp),intent(in) :: f4of2_sla(ntypat),f6of2_sla(ntypat) 137 type(pawrad_type),intent(inout) :: pawrad(ntypat) 138 type(pawtab_type),target,intent(inout) :: pawtab(ntypat) 139 140 !Local variables --------------------------------------- 141 !scalars 142 integer :: icount,il,ilmn,ilmnp,isela,iselb,itemp,itypat,iu,iup,j0lmn,jl,jlmn,jlmnp,ju,jup 143 integer :: klm0x,klma,klmb,klmn,klmna,klmnb,kln,kln1,kln2,kyc,lcur,lexexch,lkyc,ll,ll1 144 integer :: lmexexch,lmkyc,lmn_size,lmn2_size,lpawu 145 integer :: m1,m11,m2,m21,m3,m31,m4,m41 146 integer :: mesh_size,int_meshsz,mkyc,sz1 147 integer :: option_interaction_, Loc_prtvol 148 logical :: compute_euijkl,compute_euij_fll 149 real(dp) :: ak,f4of2,f6of2,int1,intg,phiint_ij,phiint_ipjp,vee1,vee2 150 character(len=500) :: message 151 152 logical,optional,intent(in) :: lmagCalc 153 logical :: lmagCalc_ 154 !arrays 155 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 156 real(dp) :: euijkl_temp(3),euijkl_temp2(3),euijkl_dc(3) 157 real(dp),allocatable :: ff(:),gg(:) 158 159 ! ************************************************************************* 160 161 DBG_ENTER("COLL") 162 Loc_prtvol = 3 163 lmagCalc_ = .False. 164 if (present(lmagCalc)) then 165 if (lmagCalc .eqv. .True.) lmagCalc_ = .True. 166 Loc_prtvol = 0 167 end if 168 169 !No correlations= nothing to do 170 if(useexexch==0.and.usepawu==0.and.use_dmft==0) then 171 do itypat=1,ntypat 172 pawtab(itypat)%usepawu=0;pawtab(itypat)%useexexch=0;pawtab(itypat)%exchmix=zero 173 end do 174 return 175 end if 176 177 !PAW+U and local exact-exchange restriction 178 if(useexexch/=0.and.usepawu/=0)then 179 do itypat=1,ntypat 180 if (llpawu(itypat)/=llexexch(itypat).and.llpawu(itypat)/=-1.and.llexexch(itypat)/=-1) then 181 write(message, '(5a,i2,3a)' )& 182 & ' When PAW+U (usepawu/=0) and local exact-exchange (exexch/=0)',ch10,& 183 & ' are selected together, they must apply on the same',ch10,& 184 & ' angular momentum (lpawu/=lexexch forbidden, here for typat=',itypat,') !',ch10,& 185 & ' Action: correct your input file.' 186 ABI_ERROR(message) 187 end if 188 end do 189 end if 190 191 !Print title 192 if((abs(usepawu)>=1.and.abs(usepawu)<=4).or.useexexch/=0.and.(.not.lmagCalc_)) & 193 & write(message, '(3a)' ) ch10,ch10," ******************************************" 194 if(usepawu==1) then 195 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL" 196 else if(usepawu==2) then 197 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: AMF" 198 else if(usepawu==3) then 199 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: AMF (alternative)" 200 else if(usepawu==4) then 201 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL with no spin polarization in the xc functional" 202 else if(usepawu==-1) then 203 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL (no use of occupation matrix) - experimental" 204 else if(usepawu==-2) then 205 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: AMF (no use of occupation matrix) - experimental" 206 else if(usepawu==-4) then 207 write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL with no spin polarization in the xc functional & 208 & (no use of occupation matrix) - experimental" 209 end if 210 if(useexexch/=0) write(message, '(3a)' ) trim(message),ch10," PAW Local Exact exchange: PBE0" 211 if((abs(usepawu)>=1.and.abs(usepawu)<=4).or.useexexch/=0 .and.(.not.lmagCalc_)) then 212 if (nspinor==2) then 213 write(message, '(3a,i1)' ) trim(message),ch10," Magnetic DC : option_interaction = ",option_interaction 214 end if 215 write(message, '(3a)' ) trim(message),ch10," ******************************************" 216 end if 217 if(use_dmft==0 .and. abs(usepawu)<=4 .and.(.not.lmagCalc_)) then 218 call wrtout(ab_out,message,'COLL') 219 call wrtout(std_out, message,'COLL') 220 end if 221 !if(use_dmft>0) then 222 !write(message, '(3a)' ) ch10, " (see DMFT data in log file) " 223 !call wrtout(ab_out,message,'COLL') 224 !endif 225 option_interaction_ = option_interaction 226 if(abs(usepawu)>=10.and.nspinor==2.and.option_interaction/=1 .and.(.not.lmagCalc_)) then 227 option_interaction_ = 1 228 write(message, '(a)' ) "When usepawu>=10, option_interaction for DC is set to 1" 229 call wrtout(std_out,message,'COLL') 230 end if 231 if(usepawu<0.and.nspinor==2.and.option_interaction_==2 .and.(.not.lmagCalc_)) then 232 write(message, '(a)' ) "option_interaction=2 is not implemented for usepawu<0. Change 'usepawu' or 'optdcmagpawu' in the input." 233 ABI_ERROR(message) 234 end if 235 236 !Loop on atom types 237 do itypat=1,ntypat 238 indlmn => pawtab(itypat)%indlmn 239 lmn_size=pawtab(itypat)%lmn_size 240 lmn2_size=pawtab(itypat)%lmn2_size 241 mesh_size=pawtab(itypat)%mesh_size 242 int_meshsz=pawrad(itypat)%int_meshsz 243 lcur=-1 244 245 ! PAW+U data 246 if (usepawu/=0.or.use_dmft>0) then 247 lcur=llpawu(itypat) 248 pawtab(itypat)%lpawu=lcur 249 if(lcur/=-1) then 250 pawtab(itypat)%usepawu=usepawu 251 pawtab(itypat)%upawu=upawu(itypat) 252 pawtab(itypat)%jpawu=jpawu(itypat) 253 pawtab(itypat)%f4of2_sla=f4of2_sla(itypat) 254 pawtab(itypat)%f6of2_sla=f6of2_sla(itypat) 255 pawtab(itypat)%option_interaction_pawu=option_interaction_ 256 else 257 pawtab(itypat)%usepawu=0 258 pawtab(itypat)%upawu=zero 259 pawtab(itypat)%jpawu=zero 260 pawtab(itypat)%f4of2_sla=zero 261 pawtab(itypat)%f6of2_sla=zero 262 pawtab(itypat)%option_interaction_pawu=option_interaction_ 263 end if 264 end if 265 266 ! Local exact-echange data 267 if (useexexch/=0) then 268 lcur=llexexch(itypat) 269 pawtab(itypat)%lexexch=lcur 270 pawtab(itypat)%exchmix=exchmix 271 if(pawtab(itypat)%lexexch==-1) pawtab(itypat)%useexexch=0 272 if(pawtab(itypat)%lexexch/=-1) pawtab(itypat)%useexexch=useexexch 273 end if 274 275 ! Select only atoms with +U 276 if(lcur/=-1) then 277 278 ! Compute number of projectors for DFT+U/local exact-exchange/DFT+DMFT 279 icount=count(indlmn(1,1:lmn_size)==lcur) 280 pawtab(itypat)%nproju=icount/(2*lcur+1) 281 if(useexexch/=0.and.pawtab(itypat)%nproju>2) then 282 write(message, '(a,a,a)' )& 283 & ' Error on the number of projectors ',ch10,& 284 & ' more than 2 projectors is not allowed for local exact-exchange' 285 ABI_ERROR(message) 286 end if 287 if(pawtab(itypat)%nproju*(2*lcur+1)/=icount) then 288 message = 'pawpuxinit: Error on the number of projectors ' 289 ABI_BUG(message) 290 end if 291 if ((.not.lmagCalc_)) then 292 write(message, '(a,a,i4,a,a,i4)' ) ch10,& 293 & ' pawpuxinit : for species ',itypat,ch10,& 294 & ' number of projectors is',pawtab(itypat)%nproju 295 call wrtout(std_out,message,'COLL') 296 297 end if 298 pawtab(itypat)%ij_proj=pawtab(itypat)%nproju*(pawtab(itypat)%nproju+1)/2 299 300 ! ================================================== 301 ! A-define useful indexes 302 ! -------------------------------------------------- 303 if (allocated(pawtab(itypat)%lnproju)) then 304 ABI_FREE(pawtab(itypat)%lnproju) 305 end if 306 ABI_MALLOC(pawtab(itypat)%lnproju,(pawtab(itypat)%nproju)) 307 icount=0 308 do ilmn=1,lmn_size 309 if(indlmn(1,ilmn)==lcur) then 310 itemp=icount/(2*lcur+1) 311 if (itemp*(2*lcur+1)==icount) then 312 pawtab(itypat)%lnproju(itemp+1)=indlmn(5,ilmn) 313 end if 314 icount=icount+1 315 end if 316 end do 317 318 if (allocated(pawtab(itypat)%klmntomn)) then 319 ABI_FREE(pawtab(itypat)%klmntomn) 320 end if 321 ABI_MALLOC(pawtab(itypat)%klmntomn,(4,lmn2_size)) 322 do jlmn=1,lmn_size 323 jl= indlmn(1,jlmn) 324 j0lmn=jlmn*(jlmn-1)/2 325 do ilmn=1,jlmn 326 il= indlmn(1,ilmn) 327 klmn=j0lmn+ilmn 328 pawtab(itypat)%klmntomn(1,klmn)=indlmn(2,ilmn)+il+1 329 pawtab(itypat)%klmntomn(2,klmn)=indlmn(2,jlmn)+jl+1 330 pawtab(itypat)%klmntomn(3,klmn)=indlmn(3,ilmn) 331 pawtab(itypat)%klmntomn(4,klmn)=indlmn(3,jlmn) 332 end do 333 end do 334 335 ! ================================================== 336 ! B-PAW+U: overlap between atomic wavefunctions 337 ! -------------------------------------------------- 338 if(dmatpuopt==1 .and.(.not.lmagCalc_)) then 339 write(message, '(4a)' ) ch10,& 340 & ' pawpuxinit : dmatpuopt=1 ',ch10,& 341 & ' PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)' 342 call wrtout(std_out,message,'COLL') 343 write(message, '(8a)' ) ch10,& 344 & ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, & 345 & ' - Is an atomic eigenfunction ',ch10, & 346 & ' - Is normalized ',ch10, & 347 & ' In other cases, choose dmatpuopt=2' 348 call wrtout(std_out,message,'COLL') 349 else if(dmatpuopt==2 .and.(.not.lmagCalc_)) then 350 write(message, '(6a)' ) ch10,& 351 & ' pawpuxinit : dmatpuopt=2 ',ch10,& 352 & ' PAW+U: dens. mat. constructed by selecting contribution',ch10,& 353 & ' for each angular momentum to the density (inside PAW augm. region(s))' 354 call wrtout(std_out,message,'COLL') 355 else if(dmatpuopt==3 .and.(.not.lmagCalc_)) then 356 write(message, '(a,a,a,a,a,a)' ) ch10,& 357 & ' pawpuxinit : dmatpuopt=3 ',ch10,& 358 & ' PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)',ch10,& 359 & ' and normalized inside PAW augm. region(s)' 360 call wrtout(std_out,message,'COLL') 361 write(message, '(6a)' ) ch10,& 362 & ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, & 363 & ' is an atomic eigenfunction',ch10, & 364 & ' In the other case, choose dmatpuopt=2' 365 call wrtout(std_out,message,'COLL') 366 end if 367 368 ABI_MALLOC(ff,(mesh_size)) 369 ff(:)=zero 370 371 if (allocated(pawtab(itypat)%ph0phiint)) then 372 ABI_FREE(pawtab(itypat)%ph0phiint) 373 end if 374 if (allocated(pawtab(itypat)%zioneff)) then 375 ABI_FREE(pawtab(itypat)%zioneff) 376 end if 377 ABI_MALLOC(pawtab(itypat)%ph0phiint,(pawtab(itypat)%nproju)) 378 ABI_MALLOC(pawtab(itypat)%zioneff,(pawtab(itypat)%nproju)) 379 380 icount=0 381 do iu=1,pawtab(itypat)%nproju 382 ! write(std_out,*)'DJA iu',iu,' mesh_size',pawtab(itypat)%mesh_size 383 ! do ju=2,pawtab(itypat)%mesh_size 384 ! ff(ju)=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawrad(itypat)%rad(ju) 385 ! write(std_out,fmt='(i5,3e15.5)')ju,pawrad(itypat)%rad(ju),ff(ju),& 386 ! & RadFnH(pawrad(itypat)%rad(ju),4,3,15.0_dp) 387 ! end do 388 ! ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))**2 389 ! call simp_gen(int1,ff,pawrad(itypat)) 390 ! write(std_out,*)'DJA iu',iu,'int1 ',int1 391 ! write(std_out,*)'DJA int1,IRadFnH',int1,IRadFnH(0.0_dp,pawrad(itypat)%rmax,4,3,12) 392 ! Calculation of zioneff 393 ju=pawtab(itypat)%mesh_size-1 394 ak=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawtab(itypat)%phi(ju+1,pawtab(itypat)%lnproju(iu)) 395 ak=ak*(pawrad(itypat)%rad(ju+1)/pawrad(itypat)%rad(ju))**(pawtab(itypat)%lpawu-1) 396 pawtab(itypat)%zioneff(iu)=log(ak)/(pawrad(itypat)%rad(ju+1)-pawrad(itypat)%rad(ju)) 397 ! Calculation of ph0phiint 398 ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(1))& 399 & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu)) 400 call simp_gen(int1,ff,pawrad(itypat)) 401 pawtab(itypat)%ph0phiint(iu)=int1 402 end do 403 if(abs(pawprtvol)>=2) then 404 do icount=1,pawtab(itypat)%nproju 405 write(message, '(a,a,i2,f9.5,a)' ) ch10,& 406 & ' pawpuxinit: icount, ph0phiint(icount)=',icount,pawtab(itypat)%ph0phiint(icount) 407 call wrtout(std_out,message,'COLL') 408 write(message, '(a,f15.5)' ) & 409 & ' pawpuxinit: zioneff=',pawtab(itypat)%zioneff(icount) 410 call wrtout(std_out,message,'COLL') 411 end do 412 write(message, '(a)' ) ch10 413 call wrtout(std_out,message,'COLL') 414 end if 415 416 if (allocated(pawtab(itypat)%phiphjint)) then 417 ABI_FREE(pawtab(itypat)%phiphjint) 418 end if 419 ABI_MALLOC(pawtab(itypat)%phiphjint,(pawtab(itypat)%ij_proj)) 420 421 icount=0 422 do ju=1,pawtab(itypat)%nproju 423 do iu=1,ju 424 icount=icount+1 425 if ((dmatpuopt==1).and.(useexexch==0)) then 426 pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)*& 427 & pawtab(itypat)%ph0phiint(ju) 428 else if((dmatpuopt==2).or.(useexexch/=0)) then 429 ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))& 430 & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(ju)) 431 call simp_gen(int1,ff,pawrad(itypat)) 432 pawtab(itypat)%phiphjint(icount)=int1 433 else if((dmatpuopt>=3).and.(useexexch==0)) then 434 pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)* & 435 & pawtab(itypat)%ph0phiint(ju)/pawtab(itypat)%ph0phiint(1)**(dmatpuopt-2) 436 else 437 write(message, '(3a)' )& 438 & ' PAW+U: dmatpuopt has a wrong value !',ch10,& 439 & ' Action : change value in input file' 440 ABI_ERROR(message) 441 end if 442 end do 443 end do 444 if(pawtab(itypat)%ij_proj/=icount) then 445 message = ' Error in the loop for calculating phiphjint ' 446 ABI_ERROR(message) 447 end if 448 ABI_FREE(ff) 449 if(abs(pawprtvol)>=2) then 450 do icount=1,pawtab(itypat)%ij_proj 451 write(message, '(a,a,i2,f9.5,a)' ) ch10,& 452 & ' PAW+U: icount, phiphjint(icount)=',icount,pawtab(itypat)%phiphjint(icount) 453 call wrtout(std_out,message,'COLL') 454 end do 455 end if 456 ! end if 457 458 ! ====================================================================== 459 ! C-PAW+U: Matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]] 460 ! 1. angular part computed from Gaunt coefficients 461 ! -------------------------------------------------------------------- 462 if (usepawu/=0) then 463 lpawu=lcur 464 465 if (allocated(pawtab(itypat)%vee)) then 466 ABI_FREE(pawtab(itypat)%vee) 467 end if 468 sz1=2*lpawu+1 469 ABI_MALLOC(pawtab(itypat)%vee,(sz1,sz1,sz1,sz1)) 470 call calc_vee(pawtab(itypat)%f4of2_sla,pawtab(itypat)%f6of2_sla,pawtab(itypat)%jpawu,& 471 & pawtab(itypat)%lpawu,pawang,pawtab(itypat)%upawu,pawtab(itypat)%vee,Loc_prtvol) 472 473 ! testu=0 474 ! write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)" 475 ! do m1=1,2*lpawu+1 476 ! write(std_out,'(2x,14(f12.6,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2),m2=1,2*lpawu+1) 477 ! do m2=1,2*lpawu+1 478 ! testu=testu+ pawtab(itypat)%vee(m1,m2,m1,m2) 479 ! enddo 480 ! enddo 481 if (.not.lmagCalc_) then 482 write(message,'(a)') ch10 483 call wrtout(std_out,message,'COLL') 484 write(message,'(a)') " Matrix of interaction vee(m1,m2,m1,m2)" 485 call wrtout(std_out,message,'COLL') 486 do m1=1,2*lpawu+1 487 write(message,'(2x,14(f20.14,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2)*Ha_eV,m2=1,2*lpawu+1) 488 call wrtout(std_out,message,'COLL') 489 ! do m2=1,2*lpawu+1 490 ! testu=testu+ pawtab(itypat)%vee(m1,m2,m1,m2) 491 ! enddo 492 enddo 493 write(message,'(a)') ch10 494 call wrtout(std_out,message,'COLL') 495 end if 496 497 ! testu=testu/((two*lpawu+one)**2) 498 ! write(std_out,*) "------------------------" 499 ! write(std_out,'(a,f12.6)') " U=", testu 500 ! write(std_out,*) "------------------------" 501 ! write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)" 502 ! do m1=1,2*lpawu+1 503 ! write(std_out,'(2x,14(f12.6,2x))') ((pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1)),m2=1,2*lpawu+1) 504 ! do m2=1,2*lpawu+1 505 ! if(m1/=m2) testumj=testumj+ pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1) 506 ! enddo 507 ! enddo 508 ! testumj=testumj/((two*lpawu)*(two*lpawu+one)) 509 ! write(std_out,*) "------------------------" 510 ! write(std_out,'(a,f12.6)') " U-J=", testumj 511 ! write(std_out,*) "------------------------" 512 ! write(std_out,*) "------------------------" 513 ! write(std_out,'(a,f12.6)') " J=", testu-testumj 514 ! write(std_out,*) "------------------------" 515 516 ! c. For DFPT (or with exp. values usepawu=-1,-2 or -4), compute euijkl 517 ! --------------------------------------------- 518 compute_euijkl=(is_dfpt.or.usepawu<0) 519 if (compute_euijkl) then 520 if (allocated(pawtab(itypat)%euijkl)) then 521 ABI_FREE(pawtab(itypat)%euijkl) 522 end if 523 ABI_MALLOC(pawtab(itypat)%euijkl,(3,lmn_size,lmn_size,lmn_size,lmn_size)) 524 pawtab(itypat)%euijkl = zero 525 compute_euij_fll = .false. 526 euijkl_temp2=zero 527 if (abs(usepawu)==1.or.abs(usepawu)==4) then ! Only for FLL 528 if (allocated(pawtab(itypat)%euij_fll)) then ! allocate euij_fll for FLL 529 ABI_FREE(pawtab(itypat)%euij_fll) 530 end if 531 ABI_MALLOC(pawtab(itypat)%euij_fll,(lmn2_size)) 532 pawtab(itypat)%euij_fll = zero 533 compute_euij_fll = .true. 534 end if 535 536 ! loop on i,j 537 do klmna=1,lmn2_size 538 ilmn=pawtab(itypat)%indklmn(7,klmna) ! i 539 jlmn=pawtab(itypat)%indklmn(8,klmna) ! j 540 if (pawtab(itypat)%indlmn(1,ilmn)==lpawu.and.pawtab(itypat)%indlmn(1,jlmn)==lpawu) then ! only correlated orbitals 541 iu = pawtab(itypat)%indlmn(3,ilmn) ! ni 542 ju = pawtab(itypat)%indlmn(3,jlmn) ! nj 543 phiint_ij = pawtab(itypat)%phiphjint(iu+(ju*(ju-1))/2) ! iu <= ju by construction (ilmn<=jlmn) 544 m2 = pawtab(itypat)%indlmn(2,ilmn) ! mi 545 m21=m2+lpawu+1 546 m1 = pawtab(itypat)%indlmn(2,jlmn) ! mj 547 m11=m1+lpawu+1 548 549 if (compute_euij_fll.and.m1==m2) then ! FLL 550 pawtab(itypat)%euij_fll(klmna) = - half * phiint_ij * ( pawtab(itypat)%jpawu - pawtab(itypat)%upawu ) 551 end if 552 553 ! loop on ip,jp (=k,l) 554 do klmnb=1,lmn2_size 555 ilmnp=pawtab(itypat)%indklmn(7,klmnb) ! ip (=k) 556 jlmnp=pawtab(itypat)%indklmn(8,klmnb) ! jp (=l) 557 if (pawtab(itypat)%indlmn(1,ilmnp)==lpawu.and.pawtab(itypat)%indlmn(1,jlmnp)==lpawu) then ! correlated orbitals 558 iup = pawtab(itypat)%indlmn(3,ilmnp) ! nip 559 jup = pawtab(itypat)%indlmn(3,jlmnp) ! njp 560 phiint_ipjp = pawtab(itypat)%phiphjint(iup+(jup*(jup-1))/2) ! iup <= jup by construction (ilmnp<=jlmnp) 561 m4 = pawtab(itypat)%indlmn(2,ilmnp) ! mip 562 m41=m4+lpawu+1 563 m3 = pawtab(itypat)%indlmn(2,jlmnp) ! mjp 564 m31=m3+lpawu+1 565 566 euijkl_dc(:) = zero 567 ! Compute the double-counting part of euijkl (invariant when exchanging i<-->j or ip<-->jp) 568 ! Must be consistent with pawuenergy and pawpupot 569 if (m1==m2.and.m3==m4) then ! In that case, we have to add the double-counting term 570 571 if (abs(usepawu)==1.and.nspinor==1) then ! FLL 572 573 euijkl_dc(1) = & 574 & phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu ) 575 euijkl_dc(2) = & 576 & phiint_ij * phiint_ipjp * pawtab(itypat)%upawu 577 578 else if (abs(usepawu)==2.and.nspinor==1) then ! AMF 579 580 euijkl_dc(1) = & 581 & two*lpawu/(two*lpawu+one) * phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu ) 582 euijkl_dc(2) = & 583 & phiint_ij * phiint_ipjp * pawtab(itypat)%upawu 584 585 else if (abs(usepawu)==4.or.nspinor>1) then ! FLL without polarization in XC or nspinor>1 586 587 euijkl_dc(1:2) = & 588 & phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - half*pawtab(itypat)%jpawu ) 589 590 ! Add term taking into account global magnetization 591 if (abs(usepawu)/=4.and.pawtab(itypat)%option_interaction_pawu==3) then 592 593 euijkl_dc(1) = euijkl_dc(1) - half * phiint_ij * phiint_ipjp * pawtab(itypat)%jpawu 594 euijkl_dc(2) = euijkl_dc(2) + half * phiint_ij * phiint_ipjp * pawtab(itypat)%jpawu 595 euijkl_dc(3) = eUijkl_dc(3) - phiint_ij * phiint_ipjp * pawtab(itypat)%jpawu 596 597 end if 598 599 end if 600 601 end if ! double-counting term 602 603 ! Array of size 3: 604 ! 1st element : coupled with up/up or down/down terms 605 ! 2nd element : coupled with up/down or down/up terms (collinear part) 606 ! 3rd element : coupled with up/down or down/up terms (non-collinear part) 607 euijkl_temp(:) = zero 608 euijkl_temp2(:) = zero 609 610 vee1 = pawtab(itypat)%vee(m11,m31,m21,m41) 611 ! Note : vee(13|24) = vee(23|14) ( so : i <--> j ) 612 ! vee(13|24) = vee(14|23) ( so : ip <--> jp ) 613 ! vee(13|24) = vee(24|13) ( so : i,ip <--> j,jp ) 614 ! Also : vee(13|24) = vee(31|42) ( so : i,j <--> ip,jp ) 615 ! ==> vee1 is invariant with respect to the permutations i <--> j , ip <--> jp and i,ip <--> j,jp 616 ! ( The term 'phiint_ij * phiint_ipjp' has the same properties) 617 euijkl_temp(1:2) = phiint_ij * phiint_ipjp * vee1 618 619 vee2 = pawtab(itypat)%vee(m11,m31,m41,m21) 620 ! Note : vee(13|42) = vee(43|12) ( so : ip <--> j ) 621 ! vee(13|42) = vee(12|43) ( so : i <--> jp ) 622 ! vee(13|42) = vee(42|13) ( so : i,ip <--> jp,j ) 623 ! Also : vee(13|42) = vee(31|24) ( so : i,j <--> ip,jp ) 624 ! Combining the third and fourth rule we get: 625 ! vee(13|42) = vee(42|13) = vee(24|31) ( so : i,ip <--> j,jp ) 626 ! ==> vee2 is invariant only with respect to the permutation i,ip <--> j,jp 627 628 ! Terms i,j,ip,jp (m2,m1,m4,m3) and j,i,jp,ip (m1,m2,m3,m4) 629 euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2 630 euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2 631 pawtab(itypat)%euijkl(:,ilmn,jlmn,ilmnp,jlmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:) 632 pawtab(itypat)%euijkl(:,jlmn,ilmn,jlmnp,ilmnp) = pawtab(itypat)%euijkl(:,ilmn,jlmn,ilmnp,jlmnp) 633 634 ! Term j,i,ip,jp (m1,m2,m4,m3) 635 vee2 = pawtab(itypat)%vee(m21,m31,m41,m11) 636 euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2 637 euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2 638 pawtab(itypat)%euijkl(:,jlmn,ilmn,ilmnp,jlmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:) 639 640 ! Term i,j,jp,ip (m2,m1,m3,m4) 641 vee2 = pawtab(itypat)%vee(m11,m41,m31,m21) 642 euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2 643 euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2 644 pawtab(itypat)%euijkl(:,ilmn,jlmn,jlmnp,ilmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:) 645 646 end if ! correlated orbitals 647 end do ! klmnb 648 end if ! correlated orbitals 649 end do ! klmna 650 651 end if ! compute_euijkl 652 end if ! usepawu 653 654 ! ====================================================================== 655 ! D-Local ex-exchange: Matrix elements of coulomb interaction and Fk 656 ! ---------------------------------------------------------------------- 657 if (useexexch/=0) then 658 lexexch=lcur 659 660 ! a. compute F(k) 661 ! --------------------------------------------- 662 if (allocated(pawtab(itypat)%fk)) then 663 ABI_FREE(pawtab(itypat)%fk) 664 end if 665 ABI_MALLOC(pawtab(itypat)%fk,(6,4)) 666 pawtab(itypat)%fk=zero 667 ABI_MALLOC(ff,(mesh_size)) 668 ABI_MALLOC(gg,(mesh_size)) 669 ff(:)=zero;gg(:)=zero 670 kln=(pawtab(itypat)%lnproju(1)*( pawtab(itypat)%lnproju(1)+1)/2) 671 do ll=1,lexexch+1 672 ll1=2*ll-2 673 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 674 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln) 675 call poisson(ff,ll1,pawrad(itypat),gg) 676 ff(1)=zero 677 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln)*gg(2:mesh_size))& 678 & /pawrad(itypat)%rad(2:mesh_size) 679 call simp_gen(intg,ff,pawrad(itypat)) 680 pawtab(itypat)%fk(1,ll)=intg*(two*ll1+one) 681 end do 682 if (pawtab(itypat)%nproju==2) then 683 kln1=kln+pawtab(itypat)%lnproju(1) 684 kln2=kln1+1 685 do ll=1,lexexch+1 686 ll1=2*ll-2 687 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 688 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1) 689 call poisson(ff,ll1,pawrad(itypat),gg) 690 ff(1)=zero 691 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))& 692 & /pawrad(itypat)%rad(2:mesh_size) 693 call simp_gen(intg,ff,pawrad(itypat)) 694 pawtab(itypat)%fk(2,ll)=intg*(two*ll1+one) 695 end do 696 do ll=1,lexexch+1 697 ll1=2*ll-2 698 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 699 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln2) 700 call poisson(ff,ll1,pawrad(itypat),gg) 701 ff(1)=zero 702 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))& 703 & /pawrad(itypat)%rad(2:mesh_size) 704 call simp_gen(intg,ff,pawrad(itypat)) 705 pawtab(itypat)%fk(3,ll)=intg*(two*ll1+one) 706 end do 707 do ll=1,lexexch+1 708 ll1=2*ll-2 709 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 710 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln) 711 call poisson(ff,ll1,pawrad(itypat),gg) 712 ff(1)=zero 713 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))& 714 & /pawrad(itypat)%rad(2:mesh_size) 715 call simp_gen(intg,ff,pawrad(itypat)) 716 pawtab(itypat)%fk(4,ll)=intg*(two*ll1+one) 717 end do 718 do ll=1,lexexch+1 719 ll1=2*ll-2 720 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 721 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln) 722 call poisson(ff,ll1,pawrad(itypat),gg) 723 ff(1)=zero 724 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))& 725 & /pawrad(itypat)%rad(2:mesh_size) 726 call simp_gen(intg,ff,pawrad(itypat)) 727 pawtab(itypat)%fk(5,ll)=intg*(two*ll1+one) 728 end do 729 do ll=1,lexexch+1 730 ll1=2*ll-2 731 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 732 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1) 733 call poisson(ff,ll1,pawrad(itypat),gg) 734 ff(1)=zero 735 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))& 736 & /pawrad(itypat)%rad(2:mesh_size) 737 call simp_gen(intg,ff,pawrad(itypat)) 738 pawtab(itypat)%fk(6,ll)=intg*(two*ll1+one) 739 end do 740 f4of2=0.6681_dp 741 f6of2=0.4943_dp 742 end if 743 ABI_FREE(ff) 744 ABI_FREE(gg) 745 746 ! b. Compute vex. 747 ! --------------------------------------------- 748 if (allocated(pawtab(itypat)%vex)) then 749 ABI_FREE(pawtab(itypat)%vex) 750 end if 751 sz1=2*lexexch+1 752 ABI_MALLOC(pawtab(itypat)%vex,(sz1,sz1,sz1,sz1,4)) 753 pawtab(itypat)%vex=zero 754 lmexexch=(lexexch-1)**2+2*(lexexch-1)+1 ! number of m value below correlated orbitals 755 klm0x=lmexexch*(lmexexch+1)/2 ! value of klmn just below correlated orbitals 756 ! --------- 4 loops for interaction matrix 757 do m1=-lexexch,lexexch 758 m11=m1+lexexch+1 759 do m2=-lexexch,m1 760 m21=m2+lexexch+1 761 ! klma= number of pair before correlated orbitals + 762 ! number of pair for m1 lower than correlated orbitals 763 ! (m1+lexexch+1)*(lexexch-1) + number of pairs for correlated orbitals 764 ! before (m1,m2) + number of pair for m2 lower than current value 765 klma=klm0x+m11*lmexexch+(m11-1)*m11/2+m21 766 do m3=-lexexch,lexexch 767 m31=m3+lexexch+1 768 do m4=-lexexch,m3 769 m41=m4+lexexch+1 770 klmb=klm0x+m31*lmexexch+(m31-1)*m31/2+m41 771 ! --------- loop on k=1,2,3 (4 if f orbitals) 772 do kyc=1,2*lexexch+1,2 773 lkyc=kyc-1 774 ll=(kyc+1)/2 775 lmkyc=(lkyc+1)*(lkyc)+1 776 ak=zero 777 do mkyc=-lkyc,lkyc,1 778 isela=pawang%gntselect(lmkyc+mkyc,klma) 779 iselb=pawang%gntselect(lmkyc+mkyc,klmb) 780 if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb) 781 end do 782 ! ----- end loop on k=1,2,3 (4 if f orbitals) 783 pawtab(itypat)%vex(m11,m31,m21,m41,ll)=ak/(two*dble(lkyc)+one) 784 end do !kyc 785 do ll=1,4 786 pawtab(itypat)%vex(m11,m31,m21,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)*four_pi 787 pawtab(itypat)%vex(m21,m31,m11,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll) 788 pawtab(itypat)%vex(m11,m41,m21,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll) 789 pawtab(itypat)%vex(m21,m41,m11,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll) 790 end do 791 end do 792 end do 793 end do 794 end do 795 796 end if !useexexch/=0 797 798 if (present(ucrpa)) then 799 if (ucrpa>=1) then 800 call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat)) 801 call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat),pawtab(itypat)%rpaw) 802 end if 803 end if 804 end if !lcur/=-1 805 end do !end loop on typat 806 807 DBG_EXIT("COLL") 808 809 end subroutine pawpuxinit
m_paw_correlations/pawuenergy [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
pawuenergy
FUNCTION
Compute contributions to energy for PAW+U calculations
INPUTS
iatom=index of current atom (absolute index, the index on current proc) noccmmp(2*lpawu+1,2*lpawu+1,nspden)=density matrix in the PAW augm. region nocctot(nspden)=number of electrons in the correlated subspace pawprtvol=control print volume and debugging output for PAW pawtab <type(pawtab_type)>=paw tabulated starting data: %lpawu=l used for dft+u %vee(2*lpawu+1*4)=screened coulomb matrix dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft= optional arguments for DMFT
OUTPUT
edftumdc= PAW+U contribution to total energy edftumdcdc= PAW+U contribution to double-counting total energy
SOURCE
1006 subroutine pawuenergy(iatom,edftumdc,edftumdcdc,noccmmp,nocctot,pawprtvol,pawtab,& 1007 & dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft) ! optional arguments (DMFT) 1008 1009 !Arguments --------------------------------------------- 1010 !scalars 1011 integer,intent(in) :: iatom,pawprtvol 1012 integer,optional,intent(in) :: dmft_dc 1013 real(dp),intent(in) :: noccmmp(:,:,:,:),nocctot(:) 1014 real(dp),intent(inout) :: edftumdc,edftumdcdc 1015 real(dp),optional,intent(inout) :: e_ee,e_dc,e_dcdc 1016 real(dp),optional,intent(in) :: j_dmft,u_dmft 1017 type(pawtab_type),intent(in) :: pawtab 1018 1019 !Local variables --------------------------------------- 1020 !scalars 1021 integer :: cplex_occ,dmftdc,ispden,jspden,lpawu,m1,m11,m2,m21,m3,m31,m4,m41,nspden 1022 real(dp) :: eks_opt3,edcdc_opt3,edcdctemp,edctemp,edftutemp,jpawu,jpawu_dc,mnorm,mx,my,mz 1023 real(dp) :: n_sig,n_sigs,n_msig,n_msigs,n_dndn,n_tot,n_upup 1024 real(dp) :: n12_ud_im,n12_du_im 1025 real(dp) :: n12_ud_re,n12_du_re 1026 real(dp) :: n34_ud_im,n34_du_im 1027 real(dp) :: n34_ud_re,n34_du_re 1028 real(dp) :: upawu 1029 real(dp),allocatable :: n12_sig(:),n34_msig(:),n34_sig(:) 1030 character(len=500) :: message 1031 1032 ! ***************************************************** 1033 1034 nspden=size(nocctot) 1035 cplex_occ=size(noccmmp,1) 1036 1037 if (size(noccmmp,4)/=nspden) then 1038 message='size of nocctot and noccmmp are inconsistent!' 1039 ABI_BUG(message) 1040 end if 1041 if (pawtab%usepawu<0) then 1042 message='not allowed for usepawu<0!' 1043 ABI_BUG(message) 1044 end if 1045 if(present(dmft_dc)) then 1046 dmftdc=dmft_dc 1047 if(pawtab%usepawu<10) then 1048 write(message,'(a,i5)') "usepawu should be =10 if dmft_dc is present ",pawtab%usepawu 1049 ABI_BUG(message) 1050 end if 1051 else 1052 dmftdc=0 1053 end if 1054 1055 DBG_ENTER("COLL") 1056 1057 lpawu=pawtab%lpawu 1058 upawu=pawtab%upawu;if(present(u_dmft)) upawu=u_dmft 1059 jpawu=pawtab%jpawu;if(present(j_dmft)) jpawu=j_dmft 1060 1061 !====================================================== 1062 !Compute DFT+U Energy 1063 !----------------------------------------------------- 1064 1065 edftutemp=zero 1066 edcdc_opt3=zero 1067 eks_opt3=zero 1068 1069 ABI_MALLOC(n12_sig,(cplex_occ)) 1070 ABI_MALLOC(n34_msig,(cplex_occ)) 1071 ABI_MALLOC(n34_sig,(cplex_occ)) 1072 do ispden=1,min(nspden,2) 1073 jspden=min(nspden,2)-ispden+1 1074 1075 ! Compute n_sigs and n_msigs for pawtab%usepawu=3 1076 if (nspden<=2) then 1077 n_sig =nocctot(ispden) 1078 n_msig=nocctot(jspden) 1079 n_tot=n_sig+n_msig 1080 else 1081 n_tot=nocctot(1) 1082 mx=nocctot(2) 1083 my=nocctot(3) 1084 mz=nocctot(4) 1085 mnorm=sqrt(mx*mx+my*my+mz*mz) 1086 if (ispden==1) then 1087 ! n_sig =half*(n_tot+mnorm) 1088 ! n_msig=half*(n_tot-mnorm) 1089 n_sig =half*(n_tot+sign(mnorm,mz)) 1090 n_msig=half*(n_tot-sign(mnorm,mz)) 1091 else 1092 ! n_sig =half*(n_tot-mnorm) 1093 ! n_msig=half*(n_tot+mnorm) 1094 n_sig =half*(n_tot-sign(mnorm,mz)) 1095 n_msig=half*(n_tot+sign(mnorm,mz)) 1096 end if 1097 end if 1098 n_sigs =n_sig/(float(2*lpawu+1)) 1099 n_msigs =n_msig/(float(2*lpawu+1)) 1100 ! if(pawtab%usepawu==3) then 1101 ! write(message,fmt=12) "noccmmp11 ",ispden,noccmmp(1,1,1,ispden) 1102 ! call wrtout(std_out,message,'COLL') 1103 ! write(message,fmt=12) "noccmmp11 ",jspden,noccmmp(1,1,1,jspden) 1104 ! call wrtout(std_out,message,'COLL') 1105 ! write(message,fmt=12) "n_sig ",ispden,n_sig 1106 ! call wrtout(std_out,message,'COLL') 1107 ! write(message,fmt=12) "n_msig ",jspden,n_msig 1108 ! call wrtout(std_out,message,'COLL') 1109 ! write(message,fmt=12) "n_sigs ",ispden,n_sigs 1110 ! call wrtout(std_out,message,'COLL') 1111 ! write(message,fmt=12) "n_msigs ",jspden,n_msigs 1112 ! call wrtout(std_out,message,'COLL') 1113 ! endif 1114 ! 12 format(a,i4,e20.10) 1115 1116 ! Compute interaction energy E_{ee} 1117 do m1=-lpawu,lpawu 1118 m11=m1+lpawu+1 1119 do m2=-lpawu,lpawu 1120 m21=m2+lpawu+1 1121 n12_sig(:)=noccmmp(:,m11,m21,ispden) 1122 if(m21==m11.and.(pawtab%usepawu==3.or.dmftdc==3)) n12_sig(1)=n12_sig(1)-n_sigs 1123 do m3=-lpawu,lpawu 1124 m31=m3+lpawu+1 1125 do m4=-lpawu,lpawu 1126 m41=m4+lpawu+1 1127 n34_sig(:) =noccmmp(:,m31,m41,ispden) 1128 n34_msig(:)=noccmmp(:,m31,m41,jspden) 1129 if(m31==m41.and.(pawtab%usepawu==3.or.dmftdc==3)) then 1130 n34_sig(1)= n34_sig(1) - n_sigs 1131 n34_msig(1)= n34_msig(1) - n_msigs 1132 end if 1133 edftutemp=edftutemp & 1134 & + n12_sig(1)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) & 1135 & + n12_sig(1)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21)) 1136 if(cplex_occ==2) then 1137 edftutemp=edftutemp & 1138 & - n12_sig(2)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) & 1139 & - n12_sig(2)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21)) 1140 end if 1141 if (pawtab%usepawu==3.or.dmftdc==3) then 1142 edcdc_opt3=edcdc_opt3 & 1143 & + n_sigs*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) & 1144 & + n_sigs*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21)) 1145 eks_opt3=eks_opt3 & 1146 & + noccmmp(1,m11,m21,ispden)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) & 1147 & + noccmmp(1,m11,m21,ispden)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21)) 1148 if(cplex_occ==2) then 1149 eks_opt3=eks_opt3 & 1150 & - noccmmp(2,m11,m21,ispden)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) & 1151 & - noccmmp(2,m11,m21,ispden)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21)) 1152 end if 1153 end if 1154 end do ! m4 1155 end do ! m3 1156 end do ! m2 1157 end do ! m1 1158 1159 end do ! ispden 1160 if (nspden==1) edftutemp=two*edftutemp ! Non-magn. system: sum up and dn energies 1161 ABI_FREE(n12_sig) 1162 ABI_FREE(n34_msig) 1163 ABI_FREE(n34_sig) 1164 1165 !Non-collinear magnetism: add non-diagonal term; see (Eq 3) in PRB 72, 024458 (2005) [[cite:Shurikov2005]] 1166 if (nspden==4) then 1167 do m1=-lpawu,lpawu 1168 m11=m1+lpawu+1 1169 do m2=-lpawu,lpawu 1170 m21=m2+lpawu+1 1171 n12_ud_re=noccmmp(1,m11,m21,3) ! updn 1172 n12_ud_im=noccmmp(2,m11,m21,3) ! updn 1173 n12_du_re=noccmmp(1,m11,m21,4) ! dnup 1174 n12_du_im=noccmmp(2,m11,m21,4) ! dnup 1175 do m3=-lpawu,lpawu 1176 m31=m3+lpawu+1 1177 do m4=-lpawu,lpawu 1178 m41=m4+lpawu+1 1179 n34_ud_re=noccmmp(1,m31,m41,3) ! updn 1180 n34_ud_im=noccmmp(2,m31,m41,3) ! updn 1181 n34_du_re=noccmmp(1,m31,m41,4) ! dnup 1182 n34_du_im=noccmmp(2,m31,m41,4) ! dnup 1183 edftutemp=edftutemp-pawtab%vee(m11,m31,m41,m21) & 1184 & *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im & 1185 & +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im) 1186 if (pawtab%usepawu==3.or.dmftdc==3) then 1187 eks_opt3=eks_opt3-pawtab%vee(m11,m31,m41,m21) & 1188 & *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im & 1189 & +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im) 1190 end if 1191 end do ! m4 1192 end do ! m3 1193 end do ! m2 1194 end do ! m1 1195 end if 1196 1197 !Divide edftutemp by 2; see (Eq 1) in PRB 77, 155104 (2008) [[cite:Amadon2008a]] 1198 edftutemp=half*edftutemp 1199 1200 !if (nspden==1) then 1201 !n_tot=two*nocctot(1) 1202 !n_upup=nocctot(1) 1203 !n_dndn=nocctot(1) 1204 !else if (nspden==2) then 1205 !n_tot=nocctot(1)+nocctot(2) 1206 !n_upup=nocctot(1) 1207 !n_dndn=nocctot(2) 1208 !else if (nspden==4) then 1209 !n_tot=nocctot(1) 1210 !mx=nocctot(2) 1211 !my=nocctot(3) 1212 !mz=nocctot(4) 1213 !mnorm=sqrt(mx*mx+my*my+mz*mz) 1214 !n_upup=half*(n_tot+mnorm) 1215 !n_dndn=half*(n_tot-mnorm) 1216 !end if 1217 n_upup=n_sig 1218 n_dndn=n_msig 1219 1220 edcdctemp=zero;edctemp=zero 1221 1222 !Full localized limit 1223 if((pawtab%usepawu==1.or.pawtab%usepawu==4).or.(dmftdc==1.or.dmftdc==4.or.dmftdc==5)) then 1224 jpawu_dc=jpawu 1225 if(dmftdc==4) then 1226 jpawu_dc=zero 1227 end if 1228 edcdctemp=edcdctemp-half*upawu*n_tot**2 1229 edctemp =edctemp +half*upawu*(n_tot*(n_tot-one)) 1230 if (nspden/=4.or.pawtab%option_interaction_pawu==2) then 1231 if(dmftdc/=5.and.pawtab%usepawu/=4) then 1232 edcdctemp=edcdctemp+half*jpawu_dc*(n_upup**2+n_dndn**2) 1233 edctemp =edctemp -half*jpawu_dc*(n_upup*(n_upup-one)+n_dndn*(n_dndn-one)) 1234 else if(dmftdc==5.or.pawtab%usepawu==4) then 1235 edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2 1236 edctemp =edctemp -quarter*jpawu_dc*(n_tot*(n_tot-two)) 1237 end if 1238 else if (nspden==4.and.(pawtab%usepawu==4.or.pawtab%option_interaction_pawu==1)) then 1239 ! write(message,'(a)') " warning: option_interaction==1 for test " 1240 ! call wrtout(std_out,message,'COLL') 1241 edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2 1242 edctemp =edctemp -quarter*jpawu_dc*(n_tot*(n_tot-two)) 1243 else if (nspden==4.and.pawtab%option_interaction_pawu==3) then 1244 ! edcdctemp= \frac{J}/{4}[ N(N) + \vect{m}.\vect{m}] 1245 edcdctemp=edcdctemp+quarter*jpawu_dc*(n_tot**2 + & 1246 & mx**2+my**2+mz**2) ! +\frac{J}/{4}\vect{m}.\vect{m} 1247 ! edctemp= -\frac{J}/{4}[ N(N-2) + \vect{m}.\vect{m}] 1248 edctemp =edctemp -quarter*jpawu_dc*( & 1249 & (n_tot*(n_tot-two)) + & 1250 & mx**2+my**2+mz**2) ! -\frac{J}/{4}\vect{m}.\vect{m} 1251 end if 1252 1253 ! Around mean field 1254 else if(pawtab%usepawu==2.or.dmftdc==2) then 1255 edctemp=edctemp+upawu*(n_upup*n_dndn)& 1256 & +half*(upawu-jpawu)*(n_upup**2+n_dndn**2) & 1257 & *(dble(2*lpawu)/dble(2*lpawu+1)) 1258 edcdctemp=-edctemp 1259 else if(pawtab%usepawu==6.or.dmftdc==6) then 1260 edctemp=edctemp+upawu*(n_tot*n_tot/4_dp)& 1261 & +half*(upawu-jpawu)*(n_tot**2+n_tot**2)/4_dp & 1262 & *(dble(2*lpawu)/dble(2*lpawu+1)) 1263 edcdctemp=-edctemp 1264 else if(pawtab%usepawu==3.or.dmftdc==3) then 1265 edcdctemp=edcdc_opt3 1266 if(abs(pawprtvol)>=3) then 1267 write(message,fmt=11) "edcdc_opt3 ",edcdc_opt3 1268 call wrtout(std_out,message,'COLL') 1269 write(message,fmt=11) "eks_opt3 ",eks_opt3 1270 call wrtout(std_out,message,'COLL') 1271 write(message,fmt=11) "eks+edcdc_opt3 ",eks_opt3+edcdc_opt3 1272 call wrtout(std_out,message,'COLL') 1273 write(message,fmt=11) "(eks+edcdc_opt3)/2 ",(eks_opt3+edcdc_opt3)/2.d0 1274 call wrtout(std_out,message,'COLL') 1275 end if 1276 end if 1277 1278 edftumdc =edftumdc +edftutemp-edctemp 1279 edftumdcdc=edftumdcdc-edftutemp-edcdctemp 1280 1281 !if(pawtab%usepawu/=10.or.pawprtvol>=3) then 1282 if(abs(pawprtvol)>=3) then 1283 if(pawtab%usepawu<10) then 1284 write(message, '(5a,i4)')ch10,'======= DFT+U Energy terms (in Hartree) ====',ch10,& 1285 & ch10,' For Atom ',iatom 1286 else if (pawtab%usepawu >= 10) then 1287 write(message, '(5a,i4)')ch10,' === DFT+U Energy terms for the DMFT occupation matrix ==',ch10,& 1288 & ch10,' For Atom ',iatom 1289 end if 1290 1291 call wrtout(std_out,message,'COLL') 1292 write(message, '(a)' )" Contributions to the direct expression of energy:" 1293 call wrtout(std_out, message,'COLL') 1294 write(message,fmt=11) " Double counting correction =",edctemp 1295 call wrtout(std_out, message,'COLL') 1296 write(message,fmt=11) " Interaction energy =",edftutemp 1297 call wrtout(std_out, message,'COLL') 1298 write(message,fmt=11) " Total DFT+U Contribution =",edftutemp-edctemp 1299 call wrtout(std_out, message,'COLL') 1300 write(message, '(a)' )' ' 1301 call wrtout(std_out, message,'COLL') 1302 write(message, '(a)' )" For the ""Double-counting"" decomposition:" 1303 call wrtout(std_out, message,'COLL') 1304 write(message,fmt=11) " DFT+U Contribution =",-edftutemp-edcdctemp 1305 call wrtout(std_out, message,'COLL') 1306 11 format(a,e20.10) 1307 if(abs(pawprtvol)>=2) then 1308 write(message,fmt=11)" edcdctemp =",edcdctemp 1309 call wrtout(std_out, message,'COLL') 1310 write(message,fmt=11)" edftumdcdc for current atom =",-edftutemp-edcdctemp 1311 call wrtout(std_out, message,'COLL') 1312 write(message, '(a)' )' ' 1313 call wrtout(std_out, message,'COLL') 1314 write(message,fmt=11)" pawuenergy: -VUKS pred =",edftumdcdc-edftumdc 1315 call wrtout(std_out, message,'COLL') 1316 end if 1317 write(message, '(a)' )' ' 1318 call wrtout(std_out, message,'COLL') 1319 end if 1320 1321 !For DMFT calculation 1322 if(present(e_ee)) e_ee=e_ee+edftutemp 1323 if(present(e_dc)) e_dc=e_dc+edctemp 1324 if(present(e_dcdc)) e_dcdc=e_dcdc+edcdctemp 1325 1326 DBG_EXIT("COLL") 1327 1328 end subroutine pawuenergy
m_paw_correlations/pawxenergy [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
pawxenergy
FUNCTION
Compute contributions to energy for PAW+ local exact exchange calculations
INPUTS
pawprtvol=control print volume and debugging output for PAW pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab <type(pawtab_type)>=paw tabulated starting data: %lexexch=l used for local exact-exchange %vex(2*lexexch+1*4)=screened coulomb matrix
SIDE EFFECTS
eexex=energy is updated with the contribution of the cuyrrent atom
SOURCE
1352 subroutine pawxenergy(eexex,pawprtvol,pawrhoij,pawtab) 1353 1354 !Arguments --------------------------------------------- 1355 !scalars 1356 integer,intent(in) :: pawprtvol 1357 real(dp),intent(inout) :: eexex 1358 type(pawrhoij_type),intent(in) :: pawrhoij 1359 type(pawtab_type),intent(in) :: pawtab 1360 1361 !Local variables --------------------------------------- 1362 !scalars 1363 integer :: irhoij,irhoij1,ispden,jrhoij,jrhoij1,klmn,klmn1,lexexch,ll,m11,m21,m31,m41,n1 1364 integer :: n2,n3,n4,nk,nn1,nn2 1365 real(dp) :: eexextemp 1366 character(len=500) :: message 1367 !arrays 1368 integer :: indn(3,3) 1369 real(dp) :: factnk(6) 1370 1371 ! ***************************************************** 1372 1373 DBG_ENTER("COLL") 1374 1375 if (pawrhoij%qphase==2) then 1376 message='pawxenergy: local exact-exchange not compatible with qphase=2!' 1377 ABI_ERROR(message) 1378 end if 1379 1380 lexexch=pawtab%lexexch 1381 if (pawtab%nproju==1) nk=1 1382 if (pawtab%nproju==2) nk=6 1383 factnk(1)=one;factnk(2)=one;factnk(3)=one 1384 factnk(4)=two;factnk(5)=two;factnk(6)=two 1385 indn(1,1)=1;indn(1,2)=4;indn(1,3)=5 1386 indn(2,1)=4;indn(2,2)=2;indn(2,3)=6 1387 indn(3,1)=5;indn(3,2)=6;indn(3,3)=3 1388 1389 !====================================================== 1390 !Compute local exact exchange Energy 1391 !----------------------------------------------------- 1392 eexextemp=zero 1393 1394 do ispden=1,pawrhoij%nspden 1395 jrhoij=1 1396 do irhoij=1,pawrhoij%nrhoijsel 1397 klmn=pawrhoij%rhoijselect(irhoij) 1398 if(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*lexexch) then 1399 m11=pawtab%klmntomn(1,klmn);m21=pawtab%klmntomn(2,klmn) 1400 n1=pawtab%klmntomn(3,klmn);n2=pawtab%klmntomn(4,klmn) 1401 nn1=(n1*n2)/2+1 1402 jrhoij1=1 1403 do irhoij1=1,pawrhoij%nrhoijsel 1404 klmn1=pawrhoij%rhoijselect(irhoij1) 1405 if(pawtab%indklmn(3,klmn1)==0.and.pawtab%indklmn(4,klmn1)==2*lexexch) then 1406 m31=pawtab%klmntomn(1,klmn1);m41=pawtab%klmntomn(2,klmn1) 1407 n3=pawtab%klmntomn(3,klmn1);n4=pawtab%klmntomn(4,klmn1) 1408 nn2=(n3*n4)/2+1 1409 do ll=1,lexexch+1 1410 eexextemp=eexextemp-pawtab%vex(m11,m31,m41,m21,ll)*pawtab%dltij(klmn)*pawtab%fk(indn(nn1,nn2),ll)& 1411 & *pawtab%dltij(klmn1)*pawrhoij%rhoijp(jrhoij,ispden)*pawrhoij%rhoijp(jrhoij1,ispden) 1412 end do 1413 end if 1414 jrhoij1=jrhoij1+pawrhoij%cplex_rhoij 1415 end do !irhoij1 1416 end if 1417 jrhoij=jrhoij+pawrhoij%cplex_rhoij 1418 end do !irhoij 1419 end do ! ispden 1420 eexextemp=eexextemp/two 1421 eexex=eexex+eexextemp*pawtab%exchmix 1422 1423 if (abs(pawprtvol)>=2) then 1424 write(message, '(a)' )" Contributions to the direct expression of energy:" 1425 call wrtout(std_out,message,'COLL') 1426 write(message,fmt='(a,e20.10,a)') " HF exchange energy =",eexextemp,ch10 1427 call wrtout(std_out,message,'COLL') 1428 end if 1429 1430 DBG_EXIT("COLL") 1431 1432 end subroutine pawxenergy
m_paw_correlations/setnoccmmp [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
setnoccmmp
FUNCTION
PAW+U only:
Compute density matrix nocc_{m,m_prime}
or
Impose value of density matrix using dmatpawu input array, then symetrize it.
noccmmp^{\sigma}_{m,m'}=\sum_{ni,nj}[\rho^{\sigma}_{ni,nj}*phiphjint_{ni,nj}]
INPUTS
compute_dmat= flag: if 1, nocc_{m,mp} is computed
dimdmat=first dimension of dmatpawu array
dmatpawu(dimdmat,dimdmat,nsppol*nspinor,natpawu)=input density matrix to be copied into noccmpp
dmatudiag= flag controlling the use of diagonalization:
0: no diagonalization of nocc_{m,mp}
1: diagonalized nocc_{m,mp} matrix is printed
2: dmatpawu matrix is expressed in the basis where nocc_(m,mp} is diagonal
impose_dmat= flag: if 1, nocc_{m,mp} is replaced by dmatpawu
indsym(4,nsym,natom)=indirect indexing array for atom labels
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
comm_atom=--optional-- MPI communicator over atoms
my_natom=number of atoms treated by current processor
natom=number of atoms in cell
natpawu=number of atoms on which PAW+U is applied
nspinor=number of spinorial components of the wavefunctions
nsppol=number of independant spin components
nsym=number of symmetry elements in space group
ntypat=number of atom types
paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
pawang <type(pawang_type)>=paw angular mesh and related data
pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
spinat(3,matom)=initial spin of each atom, in unit of hbar/2
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
typat(natom)=type for each atom
useexexch=1 if local-exact-exchange is activated
usepawu= /=0 if PAW+U is activated
OUTPUT
paw_ij(natom)%noccmmp(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol or ndij)=density matrix
NOTES
For non-collinear magnetism,
- nocc_{m,mp} is computed as:noccmmp(:,:,:,1)= n{m,mp}
noccmmp(:,:,:,2)= m_x{m,mp}
noccmmp(:,:,:,3)= m_y{m,mp}
noccmmp(:,:,:,4)= m_z{m,mp}
- but nocc_{m,mp} is stored as: noccmmp(:,:,:,1)= n^{up,up}_{m,mp}
noccmmp(:,:,:,2)= n^{dn,dn}_{m,mp}
noccmmp(:,:,:,3)= n^{up,dn}_{m,mp}
noccmmp(:,:,:,4)= n^{dn,up}_{m,mp}
We choose to have noccmmp complex when ndij=4 (ie nspinor=2)
If ndij=4 and pawspnorb=0, one could keep noccmmp real
with the n11, n22, Re(n12), Im(n21) representation, but it would
less clear to change the representation when pawspnorb is activated.
If ndij=4, nocc_{m,mp} is transformed to the Ylm basis
and then to the J, M_J basis (if cplex_dij==2)
Note that n_{m,mp}=<mp|hat(n)|m> because rhoij=<p_j|...|p_i>
SOURCE
1502 subroutine setnoccmmp(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,& 1503 & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawprtvol,pawrhoij,pawtab,& 1504 & spinat,symafm,typat,useexexch,usepawu, & 1505 & mpi_atmtab,comm_atom,l_orbmom,atom_orbmom,my_l_occmat) ! optional arguments (parallelism) and printing lorb mag 1506 1507 !Arguments --------------------------------------------- 1508 !scalars 1509 integer,intent(in) :: compute_dmat,dimdmat,dmatudiag,impose_dmat,my_natom,natom,natpawu 1510 integer,intent(in) :: nspinor,nsppol,nsym,ntypat,useexexch,usepawu 1511 integer,optional,intent(in) :: comm_atom 1512 type(pawang_type),intent(in) :: pawang 1513 integer,intent(in) :: pawprtvol 1514 !arrays 1515 integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),typat(natom) 1516 integer,optional,target,intent(in) :: mpi_atmtab(:) 1517 integer,optional,intent(in) :: l_orbmom,atom_orbmom 1518 real(dp),intent(in) :: dmatpawu(dimdmat,dimdmat,nspinor*nsppol,natpawu*impose_dmat) 1519 real(dp),intent(in) :: spinat(3,natom) 1520 type(paw_ij_type),intent(inout) :: paw_ij(my_natom) 1521 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 1522 type(pawtab_type),intent(in) :: pawtab(ntypat) 1523 1524 !Local variables --------------------------------------- 1525 !scalars 1526 integer,parameter :: limp=0 ! could become an input variable 1527 integer :: at_indx,cplex_dij,cplex_rhoij,dmatudiag_loc,iafm,iatom,iatom_tot,iatpawu,icount 1528 integer :: ilm,im1,im2,in1,in2,info,iplex,irot,ispden, irhoij,itypat,jlm,jrhoij 1529 integer :: jspden,klmn,kspden,lcur,ldim,lmax,lmin,lpawu,lwork,my_comm_atom,ndij,nmat,nspden,nsploop 1530 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used with non-collinear magnetism 1531 logical :: antiferro,my_atmtab_allocated,noccsym_error,paral_atom,use_afm 1532 ! real(dp),parameter :: invsqrt2=one/sqrt2 1533 real(dp) :: factafm,mnorm,mx,my,mz,ntot,nup,ndn,snorm,sx,sy,szm,szp 1534 character(len=4) :: wrt_mode 1535 character(len=500) :: message 1536 !arrays 1537 integer :: nsym_used(2) 1538 integer,pointer :: my_atmtab(:) 1539 real(dp) :: ro(2),sumocc(2) 1540 real(dp),allocatable :: eig(:),hdp(:,:,:),hdp2(:,:),noccmmptemp(:,:,:,:),noccmmp_tmp(:,:,:,:) 1541 real(dp),allocatable :: rwork(:),noccmmp2(:,:,:,:),nocctot2(:) 1542 complex(dpc),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:) 1543 complex(dpc),allocatable :: zhdp(:,:),zhdp2(:,:),znoccmmp_tmp(:,:),zwork(:) 1544 character(len=9),parameter :: dspin(6)= (/"up ","down ","up-up ","down-down","Re[up-dn]","Im[up-dn]"/) 1545 character(len=9),parameter :: dspinc(6)= (/"up ","down ","up-up ","down-down","up-dn ","dn-up "/) 1546 ! character(len=9),parameter :: dspinc2(6)=(/"up ","down ","dn-dn ","up-up ","dn-up ","up-dn "/) 1547 character(len=9),parameter :: dspinm(6)= (/"dn ","up i ","n ","mx ","my ","mz "/) 1548 type(coeff4_type),allocatable :: tmp_noccmmp(:) 1549 1550 real(dp),allocatable :: l_noccmmp_tmp(:,:,:,:) 1551 real(dpc),optional,allocatable :: my_l_occmat(:,:,:,:) 1552 logical :: cal_lmom 1553 integer :: atom_min,atom_max 1554 1555 !********************************************************************* 1556 1557 DBG_ENTER("COLL") 1558 !in case of calculating orbital magnetic moments, only the occupation matrix for atoms atom_orbmom and orbital l_orbmom 1559 !is calculated and returned in my_l_occmat. 1560 if (present(l_orbmom) .and. present(atom_orbmom)) then 1561 cal_lmom= .true. 1562 atom_min=atom_orbmom 1563 atom_max=atom_orbmom 1564 else 1565 cal_lmom=.false. 1566 atom_min=1 1567 atom_max=my_natom 1568 end if 1569 1570 !Tests 1571 if (my_natom>0) then 1572 if (nsppol/=paw_ij(1)%nsppol) then 1573 message='inconsistent values for nsppol!' 1574 ABI_BUG(message) 1575 end if 1576 if (compute_dmat>0) then 1577 if (pawrhoij(1)%nspden/=paw_ij(1)%nspden.and.& 1578 & pawrhoij(1)%nspden/=4.and.paw_ij(1)%nspden/=1) then 1579 message=' inconsistent values for nspden!' 1580 ABI_BUG(message) 1581 end if 1582 end if 1583 if (pawrhoij(1)%qphase==2) then 1584 message='setnoccmmp not compatible with qphase=2!' 1585 ABI_BUG(message) 1586 end if 1587 end if 1588 if (usepawu/=0.and.useexexch/=0) then 1589 message='usepawu/=0 and useexexch>0 not allowed!' 1590 ABI_BUG(message) 1591 end if 1592 if (impose_dmat/=0.and.dimdmat==0) then 1593 message='dmatpawu must be allocated when impose_dmat/=0!' 1594 ABI_BUG(message) 1595 end if 1596 if (usepawu>0.and.compute_dmat/=0.and.impose_dmat/=0.and.pawang%nsym==0) then 1597 message='pawang%zarot must be allocated!' 1598 ABI_BUG(message) 1599 end if 1600 1601 !Some inits 1602 if (usepawu==0.and.useexexch==0) return 1603 nspden=1;ndij=1;cplex_dij=1 1604 if (my_natom>0) then 1605 nspden=paw_ij(1)%nspden 1606 ndij=paw_ij(1)%ndij 1607 cplex_dij=paw_ij(1)%cplex_dij 1608 end if 1609 antiferro=(nspden==2.and.nsppol==1) 1610 use_afm=((antiferro).or.((nspden==4).and.afm_noncoll)) 1611 dmatudiag_loc=dmatudiag 1612 if (dmatudiag==2.and.(dimdmat==0.or.impose_dmat==0)) dmatudiag_loc=1 1613 1614 !Set up parallelism over atoms 1615 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 1616 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 1617 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 1618 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) !vz_d 1619 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 1620 1621 !If needed, store dmatpu in suitable format in tmp_noccmmp 1622 if (usepawu/=0.and.impose_dmat/=0) then 1623 iatpawu=0 1624 ABI_MALLOC(tmp_noccmmp,(natom)) 1625 do iatom_tot=1,natom 1626 itypat=typat(iatom_tot) 1627 lpawu=pawtab(itypat)%lpawu 1628 if (lpawu/=-1) then 1629 iatpawu=iatpawu+1 1630 if (ndij/=4) then 1631 ABI_MALLOC(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol)) 1632 tmp_noccmmp(iatom_tot)%value(1,1:2*lpawu+1,1:2*lpawu+1,1:nsppol)=& 1633 & dmatpawu(1:2*lpawu+1,1:2*lpawu+1,1:nsppol,iatpawu) 1634 else 1635 ABI_MALLOC(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,ndij)) 1636 tmp_noccmmp(iatom_tot)%value=zero 1637 if (limp==0) then ! default reading 1638 snorm=sqrt(spinat(1,iatom_tot)**2+spinat(1,iatom_tot)**2+spinat(3,iatom_tot)**2) 1639 if (snorm>tol12.and.nspden/=1) then 1640 sx=half*spinat(1,iatom_tot)/snorm 1641 sy=half*spinat(2,iatom_tot)/snorm 1642 szp=half*(one+spinat(3,iatom_tot)/snorm) 1643 szm=half*(one-spinat(3,iatom_tot)/snorm) 1644 else 1645 sx=zero;sy=zero 1646 szp=half;szm=half 1647 end if 1648 do im2=1,2*lpawu+1 1649 do im1=1,2*lpawu+1 1650 nup=dmatpawu(im1,im2,1,iatpawu);ndn=dmatpawu(im1,im2,2,iatpawu) 1651 ! if (nspden==1) tmp_noccmmp(iatom_tot)%value(1,im1,im2,1:2)=half*(nup+ndn) 1652 tmp_noccmmp(iatom_tot)%value(1,im1,im2,1)=nup*szp+ndn*szm 1653 tmp_noccmmp(iatom_tot)%value(1,im1,im2,2)=nup*szm+ndn*szp 1654 tmp_noccmmp(iatom_tot)%value(1,im1,im2,3)=(nup-ndn)*sx 1655 tmp_noccmmp(iatom_tot)%value(1,im1,im2,4)=(ndn-nup)*sy 1656 end do 1657 end do 1658 1659 else if (limp>=1) then 1660 ABI_MALLOC(noccmmp_ylm,(2*lpawu+1,2*lpawu+1,ndij)) 1661 noccmmp_ylm=czero 1662 ABI_MALLOC(noccmmp_slm,(2*lpawu+1,2*lpawu+1,ndij)) 1663 noccmmp_slm=czero 1664 ABI_MALLOC(noccmmp_jmj,(2*(2*lpawu+1),2*(2*lpawu+1))) 1665 noccmmp_jmj=czero 1666 if(limp==1) then ! read input matrix in J,M_J basis (l-1/2, then l+1/2) 1667 noccmmp_jmj=czero 1668 do im1=1,2*lpawu+1 1669 noccmmp_jmj(im1,im1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp) 1670 noccmmp_jmj(im1+lpawu,im1+lpawu)=cmplx(dmatpawu(im1+lpawu,im1+lpawu,2,iatpawu),zero,kind=dp) 1671 end do 1672 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 1673 & ' == Imposed occupation matrix (in the J M_J basis: L-1/2 and L+1/2 states)' 1674 call wrtout(std_out,message,wrt_mode) 1675 call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,& 1676 & 2,2,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 1677 end if 1678 if (limp==2) then ! read input matrix in Ylm basis 1679 noccmmp_ylm=czero 1680 do im1=1,2*lpawu+1 1681 noccmmp_ylm(im1,im1,1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp) 1682 noccmmp_ylm(im1,im1,2)=cmplx(dmatpawu(im1,im1,2,iatpawu),zero,kind=dp) 1683 end do 1684 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 1685 & ' == Imposed occupation matrix (in the Ylm basis), for dn and up spin' 1686 call wrtout(std_out,message,wrt_mode) 1687 end if 1688 call mat_slm2ylm(lpawu,noccmmp_ylm,noccmmp_slm,ndij,& 1689 & 2,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first 1690 ! interchange upup and dndn 1691 if (limp>=1) then 1692 tmp_noccmmp(iatom_tot)%value(1,:,:,1)=real(noccmmp_slm(:,:,2)) 1693 tmp_noccmmp(iatom_tot)%value(2,:,:,1)=aimag(noccmmp_slm(:,:,2)) 1694 tmp_noccmmp(iatom_tot)%value(1,:,:,2)=real(noccmmp_slm(:,:,1)) 1695 tmp_noccmmp(iatom_tot)%value(2,:,:,2)=aimag(noccmmp_slm(:,:,1)) 1696 tmp_noccmmp(iatom_tot)%value(1,:,:,3)=real(noccmmp_slm(:,:,4)) 1697 tmp_noccmmp(iatom_tot)%value(2,:,:,3)=aimag(noccmmp_slm(:,:,4)) 1698 tmp_noccmmp(iatom_tot)%value(1,:,:,4)=real(noccmmp_slm(:,:,3)) 1699 tmp_noccmmp(iatom_tot)%value(2,:,:,4)=aimag(noccmmp_slm(:,:,3)) 1700 end if 1701 if(abs(pawprtvol)>2) then 1702 write(message, '(2a)' ) ch10,& 1703 & " Check Imposed density matrix in different basis" 1704 call wrtout(std_out,message,wrt_mode) 1705 call mat_slm2ylm(lpawu,noccmmp_slm,noccmmp_ylm,ndij,& 1706 & 1,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first 1707 call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,1,2,& 1708 & pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 1709 end if 1710 ABI_FREE(noccmmp_ylm) 1711 ABI_FREE(noccmmp_jmj) 1712 ABI_FREE(noccmmp_slm) 1713 end if 1714 end if 1715 end if 1716 end do 1717 end if ! impose_dmat/=0 1718 1719 !Print message 1720 if (usepawu/=0.and.impose_dmat/=0) then 1721 if (dmatudiag_loc/=2) then 1722 write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is kept constant',ch10,& 1723 & 'and equal to dmatpawu from input file !',ch10,& 1724 & '----------------------------------------------------------' 1725 else 1726 write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is imposed',ch10,& 1727 & 'and equal to dmatpawu in the diagonal basis !',ch10,& 1728 & '----------------------------------------------------------' 1729 end if 1730 call wrtout(std_out,message,'COLL') 1731 end if 1732 1733 if (usepawu/=0.and.dmatudiag_loc/=0.and.my_natom>0) then 1734 write(message,'(4a)') ch10,'Diagonalized occupation matrix "noccmmp" is printed !',ch10,& 1735 & '-------------------------------------------------------------' 1736 call wrtout(std_out,message,wrt_mode) 1737 end if 1738 1739 !Loops over atoms 1740 do iatom=atom_min,atom_max 1741 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 1742 itypat=pawrhoij(iatom)%itypat 1743 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 1744 1745 if (.not. cal_lmom) then 1746 if (useexexch/=0) then 1747 lcur=pawtab(itypat)%lexexch 1748 else if (usepawu/=0) then 1749 lcur=pawtab(itypat)%lpawu 1750 end if 1751 end if 1752 1753 if (cal_lmom) lcur=l_orbmom 1754 if (lcur/=-1) then 1755 1756 ! ######################################################################################## 1757 ! # Compute nocc_mmp 1758 ! ######################################################################################## 1759 if ((usepawu/=0.and.compute_dmat/=0).or.useexexch/=0) then 1760 1761 1762 ABI_MALLOC(l_noccmmp_tmp,(cplex_dij,2*lcur+1,2*lcur+1,ndij)) 1763 l_noccmmp_tmp(:,:,:,:)=zero 1764 1765 ! Loop over spin components 1766 ABI_MALLOC(noccmmptemp,(cplex_dij,2*lcur+1,2*lcur+1,ndij)) 1767 noccmmptemp(:,:,:,:)=zero 1768 if(ndij==4) then 1769 ABI_MALLOC(noccmmp2,(cplex_dij,2*lcur+1,2*lcur+1,ndij)) 1770 end if 1771 if(ndij==4) then 1772 if(allocated(nocctot2)) then 1773 ABI_FREE(nocctot2) 1774 end if 1775 ABI_MALLOC(nocctot2,(ndij)) 1776 end if 1777 nsploop=ndij 1778 do ispden=1,nsploop 1779 jrhoij=1 1780 do irhoij=1,pawrhoij(iatom)%nrhoijsel 1781 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 1782 im1=pawtab(itypat)%klmntomn(1,klmn) 1783 im2=pawtab(itypat)%klmntomn(2,klmn) 1784 in1=pawtab(itypat)%klmntomn(3,klmn) 1785 in2=pawtab(itypat)%klmntomn(4,klmn) 1786 lmin=pawtab(itypat)%indklmn(3,klmn) 1787 lmax=pawtab(itypat)%indklmn(4,klmn) 1788 1789 ro(1:2)=zero 1790 ro(1:cplex_rhoij)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+cplex_rhoij-1,ispden) 1791 if (ndij==1) ro(1:2)=half*ro(1:2) 1792 ! Non-collinear magnetism: keep n, m storage because 1793 ! it is easier for the computation of noccmmp from rhoij) 1794 1795 if(lmin==0.and.lmax==2*lcur) then 1796 icount=in1+(in2*(in2-1))/2 1797 if(pawtab(itypat)%ij_proj<icount .and. (.not. cal_lmom) ) then 1798 message='PAW+U: Problem in the loop calculating noccmmp!' 1799 ABI_BUG(message) 1800 end if 1801 if(in1/=in2) then 1802 if(im2<=im1) then 1803 noccmmptemp(1:cplex_dij,im1,im2,ispden)=noccmmptemp(1:cplex_dij,im1,im2,ispden) & 1804 & +ro(1:cplex_dij)*pawtab(itypat)%phiphjint(icount) 1805 end if 1806 end if 1807 if(im2>=im1) then 1808 l_noccmmp_tmp(1:cplex_dij,im1,im2,ispden)=l_noccmmp_tmp(1:cplex_dij,im1,im2,ispden) & 1809 & +ro(1:cplex_dij)*pawtab(itypat)%phiphjint(icount) 1810 end if 1811 end if 1812 jrhoij=jrhoij+cplex_rhoij 1813 end do ! irhoij 1814 do im2=1,2*lcur+1 1815 do im1=1,im2 1816 l_noccmmp_tmp(1,im1,im2,ispden)=l_noccmmp_tmp(1,im1,im2,ispden) & 1817 & +noccmmptemp(1,im2,im1,ispden) 1818 if(cplex_dij==2) l_noccmmp_tmp(2,im1,im2,ispden)=l_noccmmp_tmp(2,im1,im2,ispden) & 1819 & -noccmmptemp(2,im2,im1,ispden) 1820 end do 1821 end do 1822 do im1=1,2*lcur+1 1823 do im2=1,im1 1824 l_noccmmp_tmp(1,im1,im2,ispden)=l_noccmmp_tmp(1,im2,im1,ispden) 1825 if(cplex_dij==2) l_noccmmp_tmp(2,im1,im2,ispden)=-l_noccmmp_tmp(2,im2,im1,ispden) 1826 end do 1827 end do 1828 end do ! ispden 1829 ABI_FREE(noccmmptemp) 1830 ! Compute noccmmp2, occupation matrix in the spin basis (upup, dndn, updn, dnup) 1831 if(ndij==4) then 1832 noccmmp2(:,:,:,:)=zero 1833 do im1=1,2*lcur+1 1834 do im2=1,2*lcur+1 1835 noccmmp2(1,im1,im2,1)=half*(l_noccmmp_tmp(1,im1,im2,1)+l_noccmmp_tmp(1,im1,im2,4)) 1836 noccmmp2(2,im1,im2,1)=half*(l_noccmmp_tmp(2,im1,im2,1)+l_noccmmp_tmp(2,im1,im2,4)) 1837 noccmmp2(1,im1,im2,2)=half*(l_noccmmp_tmp(1,im1,im2,1)-l_noccmmp_tmp(1,im1,im2,4)) 1838 noccmmp2(2,im1,im2,2)=half*(l_noccmmp_tmp(2,im1,im2,1)-l_noccmmp_tmp(2,im1,im2,4)) 1839 noccmmp2(1,im1,im2,3)=half*(l_noccmmp_tmp(1,im1,im2,2)+l_noccmmp_tmp(2,im1,im2,3)) 1840 noccmmp2(2,im1,im2,3)=half*(l_noccmmp_tmp(2,im1,im2,2)-l_noccmmp_tmp(1,im1,im2,3)) 1841 noccmmp2(1,im1,im2,4)=half*(l_noccmmp_tmp(1,im1,im2,2)-l_noccmmp_tmp(2,im1,im2,3)) 1842 noccmmp2(2,im1,im2,4)=half*(l_noccmmp_tmp(2,im1,im2,2)+l_noccmmp_tmp(1,im1,im2,3)) 1843 end do 1844 end do 1845 if(abs(pawprtvol)>=1 .and. (.not. cal_lmom)) then 1846 write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals in the n, m basis :" 1847 call wrtout(std_out,message,wrt_mode) 1848 do ispden=1,ndij 1849 write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinm(ispden+2*(ndij/4))) 1850 call wrtout(std_out,message,wrt_mode) 1851 do im1=1,lcur*2+1 ! ( order of indices in noccmmp is exchanged in order to have the same convention as rhoij: transposition is done after ) 1852 if(cplex_dij==1)& 1853 & write(message,'(12(1x,9(1x,f10.5)))')& 1854 & (l_noccmmp_tmp(1,im2,im1,ispden),im2=1,lcur*2+1) 1855 if(cplex_dij==2)& 1856 ! & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 1857 & write(message,'(12(1x,9(1x,"(",f10.5,",",f10.5,")")))')& 1858 & (l_noccmmp_tmp(:,im2,im1,ispden),im2=1,lcur*2+1) 1859 call wrtout(std_out,message,wrt_mode) 1860 end do 1861 end do 1862 end if ! pawprtvol >=1 1863 end if 1864 1865 ! Compute total number of electrons per spin 1866 if (.not. cal_lmom) then 1867 paw_ij(iatom)%nocctot(:)=zero ! contains nmmp in the n m representation 1868 if(ndij==4) nocctot2(:)=zero ! contains nmmp in the upup dndn updn dnup representation 1869 do ispden=1,ndij 1870 do im1=1,2*lcur+1 1871 if(ndij==4) then 1872 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+l_noccmmp_tmp(1,im1,im1,ispden) 1873 nocctot2(ispden)=nocctot2(ispden)+noccmmp2(1,im1,im1,ispden) 1874 else 1875 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+l_noccmmp_tmp(1,im1,im1,ispden) 1876 end if 1877 end do 1878 end do 1879 end if 1880 ! noccmmp will now be in the up up , dn dn... representation and now n_mmp=<m|n|mp> instead of <mp|n|m> ! 1881 if(ndij==4) then 1882 do ispden=1,ndij 1883 do iplex=1,cplex_dij 1884 do im1=1,2*lcur+1 1885 do im2=1,2*lcur+1 1886 l_noccmmp_tmp(iplex,im1,im2,ispden)=noccmmp2(iplex,im2,im1,ispden) ! now, noccmmp is in the upup dndn updn dnup representation 1887 end do 1888 end do 1889 end do 1890 end do 1891 ABI_FREE(noccmmp2) 1892 end if 1893 ! Printing of new nocc_mmp 1894 if (.not. cal_lmom) then 1895 if ((usepawu/=0.and.abs(usepawu)<10).or.(usepawu>=10.and.pawprtvol>=3)) then 1896 write(message, '(2a)' ) ch10, & 1897 & '========== DFT+U DATA =================================================== ' 1898 end if 1899 if (useexexch/=0) then 1900 write(message, '(2a)' ) ch10, & 1901 & '======= Local ex-exchange (PBE0) DATA =================================== ' 1902 end if 1903 if (((usepawu/=0.and.abs(usepawu)<10).or.(usepawu>=10.and.pawprtvol>=3)).or.useexexch/=0) then 1904 call wrtout(std_out,message,wrt_mode) 1905 end if 1906 if (usepawu>=10.and.pawprtvol>=3) then 1907 write(message, '(6a)' ) ch10,' ( A DFT+DMFT calculation is carried out ',& 1908 ch10,' Thus, the following DFT+U occupation matrices are not physical ',& 1909 ch10,' and just informative )' 1910 call wrtout(std_out,message,wrt_mode) 1911 end if 1912 if(abs(usepawu)<10.or.pawprtvol>=3) then ! Always write except if DMFT and pawprtvol low 1913 write(message,'(2a,i5,a,i4,a)') ch10,"====== For Atom", iatom_tot,& 1914 & ", occupations for correlated orbitals. l =",lcur,ch10 1915 call wrtout(std_out,message,wrt_mode) 1916 if(ndij==2) then 1917 do ispden=1,2 1918 write(message,'(a,i4,3a,f10.5)') "Atom", iatom_tot,". Occupations for spin ",& 1919 & trim(dspin(ispden))," =",paw_ij(iatom)%nocctot(ispden) 1920 call wrtout(std_out,message,wrt_mode) 1921 end do 1922 write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. is ",& 1923 & paw_ij(iatom)%nocctot(2)-paw_ij(iatom)%nocctot(1) 1924 call wrtout(std_out,message,wrt_mode) 1925 end if 1926 if(ndij==4) then 1927 ntot=paw_ij(iatom)%nocctot(1) 1928 mx=paw_ij(iatom)%nocctot(2) 1929 my=paw_ij(iatom)%nocctot(3) 1930 mz=paw_ij(iatom)%nocctot(4) 1931 mnorm=sqrt(mx*mx+my*my+mz*mz) 1932 nup=nocctot2(1) 1933 ndn=nocctot2(2) 1934 write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. x is ",mx 1935 call wrtout(std_out,message,wrt_mode) 1936 write(message,'(14x,a,2x,e16.8)') " local Mag. y is ",my 1937 call wrtout(std_out,message,wrt_mode) 1938 write(message,'(14x,a,2x,e16.8)') " local Mag. z is ",mz 1939 call wrtout(std_out,message,wrt_mode) 1940 write(message,'(14x,a,2x,e16.8)') " norm of Mag. is ",mnorm 1941 call wrtout(std_out,message,wrt_mode) 1942 write(message,'(14x,a,2x,f10.5)') " occ. of majority spin is ",half*(ntot+mnorm) ! to be checked versus direct calc from noccmmp 1943 call wrtout(std_out,message,wrt_mode) 1944 if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') " occ. for spin up (along z) ",nup 1945 if(abs(pawprtvol)>=1) then 1946 call wrtout(std_out,message,wrt_mode) 1947 end if 1948 write(message,'(14x,a,2x,f10.5)') " occ. of minority spin is ",half*(ntot-mnorm) 1949 call wrtout(std_out,message,wrt_mode) 1950 if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') " occ. for spin dn (along z) ",ndn 1951 if(abs(pawprtvol)>=1) then 1952 call wrtout(std_out,message,wrt_mode) 1953 end if 1954 if(ndij==4) then 1955 ABI_FREE(nocctot2) 1956 end if 1957 end if 1958 write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals:" 1959 call wrtout(std_out,message,wrt_mode) 1960 do ispden=1,ndij 1961 write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4))) 1962 call wrtout(std_out,message,wrt_mode) 1963 do im1=1,lcur*2+1 1964 if(cplex_dij==1)& 1965 & write(message,'(12(1x,9(1x,f10.5)))')& 1966 & (l_noccmmp_tmp(1,im1,im2,ispden),im2=1,lcur*2+1) 1967 if(cplex_dij==2)& 1968 & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 1969 & (l_noccmmp_tmp(:,im1,im2,ispden),im2=1,lcur*2+1) 1970 call wrtout(std_out,message,wrt_mode) 1971 end do 1972 end do 1973 end if 1974 1975 ! Transformation matrices: real->complex spherical harmonics (for test) 1976 if(ndij==4.and.abs(pawprtvol)>=0) then 1977 ABI_MALLOC(noccmmp_ylm,(2*lcur+1,2*lcur+1,ndij)) 1978 noccmmp_ylm=czero 1979 ABI_MALLOC(noccmmp_slm,(2*lcur+1,2*lcur+1,ndij)) 1980 noccmmp_slm=czero 1981 ABI_MALLOC(noccmmp_jmj,(2*(2*lcur+1),2*(2*lcur+1))) 1982 noccmmp_jmj=czero 1983 ! go from real notation for complex noccmmp to complex notation in noccmmp_slm 1984 noccmmp_slm(:,:,:)=cmplx(l_noccmmp_tmp(1,:,:,:)& 1985 & ,l_noccmmp_tmp(2,:,:,:),kind=dp) 1986 call mat_slm2ylm(lcur,noccmmp_slm,noccmmp_ylm,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 1987 1988 do ispden=1,ndij 1989 write(message,'(3a)') ch10,"Calculated Ylm occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4))) 1990 call wrtout(std_out,message,wrt_mode) 1991 do im1=1,lcur*2+1 1992 write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') (noccmmp_ylm(im1,im2,ispden),im2=1,lcur*2+1) 1993 call wrtout(std_out,message,wrt_mode) 1994 end do 1995 end do 1996 call mat_mlms2jmj(lcur,noccmmp_ylm,noccmmp_jmj,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 1997 ABI_FREE(noccmmp_ylm) 1998 ABI_FREE(noccmmp_jmj) 1999 ABI_FREE(noccmmp_slm) 2000 end if !ndij==4 2001 paw_ij(iatom)%noccmmp(:,:,:,:)=zero 2002 paw_ij(iatom)%noccmmp(:,:,:,:)=l_noccmmp_tmp(:,:,:,:) 2003 ABI_FREE(l_noccmmp_tmp) 2004 else 2005 2006 if(allocated(my_l_occmat)) then 2007 ABI_FREE(my_l_occmat) 2008 end if 2009 2010 if(allocated(nocctot2)) then 2011 ABI_FREE(nocctot2) 2012 end if 2013 ABI_MALLOC(my_l_occmat,(cplex_dij,2*lcur+1,2*lcur+1,ndij)) 2014 my_l_occmat(:,:,:,:)=zero 2015 my_l_occmat=l_noccmmp_tmp(:,:,:,:) 2016 ABI_FREE(l_noccmmp_tmp) 2017 2018 end if ! not cal_lmom 2019 2020 end if ! impose_dmat==0 2021 2022 ! ######################################################################################## 2023 ! # Diagonalize nocc_mmp 2024 ! ######################################################################################## 2025 if(usepawu/=0.and.dmatudiag_loc>0.and.(.not. cal_lmom)) then 2026 2027 lpawu=lcur;ldim=2*lpawu+1 2028 ABI_MALLOC(noccmmp_tmp,(1,ldim,ldim,ndij)) 2029 if (ndij==4) then 2030 ABI_MALLOC(znoccmmp_tmp,(2*ldim,2*ldim)) 2031 end if 2032 2033 ! Select noccmmp for this atom 2034 do ispden=1,ndij 2035 noccmmp_tmp(1,:,:,ispden)=paw_ij(iatom)%noccmmp(1,:,:,ispden) 2036 end do 2037 if (ndij==4) then 2038 do im2=1,ldim 2039 do im1=1,ldim 2040 znoccmmp_tmp(im1 , im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1)& 2041 & ,paw_ij(iatom)%noccmmp(2,im1,im2,1),kind=dp) 2042 znoccmmp_tmp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2)& 2043 & ,paw_ij(iatom)%noccmmp(2,im1,im2,2),kind=dp) 2044 znoccmmp_tmp( im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3)& 2045 & ,paw_ij(iatom)%noccmmp(2,im1,im2,3),kind=dp) 2046 znoccmmp_tmp(ldim+im1, im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,4)& 2047 & ,paw_ij(iatom)%noccmmp(2,im1,im2,4),kind=dp) 2048 end do 2049 end do 2050 end if 2051 2052 ! Diagonalize nocc_mmp 2053 if (ndij/=4) then 2054 ABI_MALLOC(hdp,(ldim,ldim,ndij)) 2055 hdp=zero 2056 lwork=3*ldim-1 2057 ABI_MALLOC(rwork,(lwork)) 2058 ABI_MALLOC(eig,(ldim)) 2059 do ispden=1,ndij 2060 call dsyev('v','u',ldim,noccmmp_tmp(1,:,:,ispden),ldim,eig,rwork,lwork,info) 2061 if(info/=0) then 2062 message=' Error in diagonalization of noccmmp (DSYEV)!' 2063 ABI_ERROR(message) 2064 end if 2065 do ilm=1,ldim 2066 hdp(ilm,ilm,ispden)=eig(ilm) 2067 end do 2068 end do ! ispden 2069 ABI_FREE(rwork) 2070 ABI_FREE(eig) 2071 else 2072 ABI_MALLOC(hdp,(2*ldim,2*ldim,1)) 2073 hdp=zero 2074 lwork=4*ldim-1 2075 ABI_MALLOC(rwork,(6*ldim-2)) 2076 ABI_MALLOC(zwork,(lwork)) 2077 ABI_MALLOC(eig,(2*ldim)) 2078 call zheev('v','u',2*ldim,znoccmmp_tmp,2*ldim,eig,zwork,lwork,rwork,info) 2079 if(info/=0) then 2080 message=' Error in diagonalization of znoccmmp_tmp (zheev) !' 2081 ABI_ERROR(message) 2082 end if 2083 do ilm=1,2*ldim 2084 hdp(ilm,ilm,1)=eig(ilm) 2085 end do 2086 ABI_FREE(rwork) 2087 ABI_FREE(zwork) 2088 ABI_FREE(eig) 2089 end if 2090 2091 ! Print diagonalized matrix and eigenvectors 2092 do ispden=1,size(hdp,3) 2093 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Diagonalized Occupation matrix' 2094 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 2095 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 2096 if (ndij==4) write(message,fmt='(2a,i3,a)')trim(message)," ==" 2097 call wrtout(std_out,message,wrt_mode) 2098 do ilm=1,size(hdp,1) 2099 write(message,'(12(1x,9(1x,f10.5)))') (hdp(ilm,jlm,ispden),jlm=1,size(hdp,2)) 2100 call wrtout(std_out,message,wrt_mode) 2101 end do 2102 end do ! ispden 2103 if(abs(pawprtvol)>=1) then 2104 if (ndij/=4) then 2105 do ispden=1,ndij 2106 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors' 2107 if (ndij==1) write(message,fmt='(2a)') trim(message),' for spin up ==' 2108 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message),' for spin ',ispden,' ==' 2109 call wrtout(std_out,message,wrt_mode) 2110 do ilm=1,ldim 2111 write(message,'(12(1x,9(1x,f10.5)))') (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim) 2112 call wrtout(std_out,message,wrt_mode) 2113 end do 2114 end do 2115 else 2116 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors (spinors) in the real harmonics basis ==' 2117 call wrtout(std_out,message,wrt_mode) 2118 do ilm=1,2*ldim 2119 write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (znoccmmp_tmp(ilm,jlm),jlm=1,2*ldim) 2120 call wrtout(std_out,message,wrt_mode) 2121 end do 2122 end if 2123 end if 2124 2125 ! Back rotation of diagonalized matrix and printing 2126 if(abs(pawprtvol)>=1) then 2127 if (ndij/=4) then 2128 ABI_MALLOC(hdp2,(ldim,ldim)) 2129 do ispden=1,ndij 2130 call dgemm('n','t',ldim,ldim,ldim,one,hdp(:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim) 2131 call dgemm('n','n',ldim,ldim,ldim,one,noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,hdp(:,:,ispden),ldim) 2132 noccmmp_tmp(1,:,:,ispden)=hdp(:,:,ispden) 2133 end do ! ispden 2134 ABI_FREE(hdp2) 2135 else 2136 ABI_MALLOC(zhdp,(2*ldim,2*ldim)) 2137 ABI_MALLOC(zhdp2,(2*ldim,2*ldim)) 2138 zhdp(:,:)=cmplx(hdp(:,:,1),zero,kind=dp) 2139 zhdp2(:,:)=cmplx(zero,zero,kind=dp) 2140 call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim) 2141 zhdp(:,:)=cmplx(zero,zero,kind=dp) 2142 call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim) 2143 znoccmmp_tmp=zhdp 2144 ABI_FREE(zhdp) 2145 ABI_FREE(zhdp2) 2146 end if 2147 nmat=ndij ; if(ndij==4.and.cplex_dij==2) nmat=1 2148 do ispden=1,nmat 2149 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 2150 & ' == Rotated back diagonalized matrix' 2151 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 2152 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 2153 if (ndij==4.and.cplex_dij==2) write(message,fmt='(4a)') trim(message)," for all component " 2154 call wrtout(std_out,message,wrt_mode) 2155 do ilm=1,ldim*cplex_dij 2156 if(ndij==1.or.ndij==2)& 2157 & write(message,'(12(1x,9(1x,f10.5)))')& 2158 & (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim) 2159 if(ndij==4.and.cplex_dij==2)& 2160 & write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')& 2161 & (znoccmmp_tmp(ilm,jlm),jlm=1,ldim*cplex_dij) 2162 call wrtout(std_out,message,wrt_mode) 2163 end do 2164 end do ! ispden 2165 end if 2166 ABI_FREE(hdp) 2167 2168 end if ! dmatudiag_loc 2169 2170 ! ######################################################################################## 2171 ! # Impose value of nocc_mmp from dmatpu; symetrize it 2172 ! ######################################################################################## 2173 if (usepawu/=0.and.impose_dmat/=0.and.(.not.cal_lmom)) then 2174 2175 lpawu=lcur 2176 nsploop=nsppol;if (ndij==4) nsploop=4 2177 noccsym_error=.false. 2178 2179 ! Loop over spin components 2180 do ispden=1,nsploop 2181 if (ndij/=4) then 2182 jspden=min(3-ispden,paw_ij(iatom)%nsppol) 2183 else if (ispden<=2) then 2184 jspden=3-ispden 2185 else 2186 jspden=ispden 2187 end if 2188 2189 ! Loops over components of nocc_mmp 2190 do jlm=1,2*lpawu+1 2191 do ilm=1,2*lpawu+1 2192 2193 if(nsym>1.and.ndij<4) then 2194 2195 nsym_used(1:2)=0 2196 sumocc(1:2)=zero 2197 2198 ! Accumulate values of nocc_mmp over symmetries 2199 do irot=1,nsym 2200 if ((symafm(irot)/=1).and.(.not.use_afm)) cycle 2201 kspden=ispden;if (symafm(irot)==-1) kspden=jspden 2202 factafm=one;if (ispden>3) factafm=dble(symafm(irot)) 2203 iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2 2204 nsym_used(iafm)=nsym_used(iafm)+1 2205 at_indx=indsym(4,irot,iatom_tot) 2206 do im2=1,2*lpawu+1 2207 do im1=1,2*lpawu+1 2208 ! Be careful: use here R_rel^-1 in term of spherical harmonics 2209 ! which is tR_rec in term of spherical harmonics 2210 ! so, use transpose[zarot] 2211 sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(1,im1,im2,kspden) & 2212 & *pawang%zarot(im1,ilm,lpawu+1,irot)& 2213 & *pawang%zarot(im2,jlm,lpawu+1,irot) 2214 ! sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(im1,im2,kspden) & 2215 ! & *pawang%zarot(ilm,im1,lpawu+1,irot)& 2216 ! & *pawang%zarot(jlm,im2,lpawu+1,irot) 2217 end do 2218 end do 2219 end do ! End loop over symmetries 2220 2221 ! Store new values of nocc_mmp 2222 paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden)=sumocc(1)/nsym_used(1) 2223 if (.not.noccsym_error)& 2224 & noccsym_error=(abs(paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden) & 2225 & -tmp_noccmmp(iatom_tot)%value(1,ilm,jlm,ispden))>tol5) 2226 2227 ! Antiferromagnetic case: has to fill up "down" component of nocc_mmp 2228 if (antiferro.and.nsym_used(2)>0) paw_ij(iatom)%noccmmp(1,ilm,jlm,2)=sumocc(2)/nsym_used(2) 2229 2230 else ! nsym=1 2231 2232 ! Case without symetries 2233 paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden)= tmp_noccmmp(iatom_tot)%value(:,ilm,jlm,ispden) 2234 end if 2235 2236 end do !ilm 2237 end do !jlm 2238 end do ! ispden 2239 do ispden=1,nsploop 2240 paw_ij(iatom)%nocctot(ispden)=zero ! contains nmmp in the n m representation 2241 do im1=1,2*lcur+1 2242 if(ndij==4.and.ispden==1) then 2243 ! in this case, on computes total number or electron for double counting correction 2244 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 2245 & paw_ij(iatom)%noccmmp(1,im1,im1,1)+paw_ij(iatom)%noccmmp(1,im1,im1,2) 2246 else if(ndij==4.and.ispden==2) then 2247 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 2248 & paw_ij(iatom)%noccmmp(1,im1,im1,3)+paw_ij(iatom)%noccmmp(1,im1,im1,4) 2249 else if(ndij==4.and.ispden==3) then 2250 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)-& 2251 & paw_ij(iatom)%noccmmp(2,im1,im1,3)+paw_ij(iatom)%noccmmp(2,im1,im1,4) 2252 else if(ndij==4.and.ispden==4) then 2253 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 2254 & paw_ij(iatom)%noccmmp(2,im1,im1,1)-paw_ij(iatom)%noccmmp(2,im1,im1,2) 2255 else 2256 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 2257 & paw_ij(iatom)%noccmmp(1,im1,im1,ispden) 2258 end if 2259 end do 2260 end do ! ispden 2261 2262 ! Printing of new nocc_mmp 2263 do ispden=1,ndij 2264 if(dmatudiag_loc==2) then 2265 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 2266 & ' == Imposed occupation matrix (in the basis of diagonalization!!)' 2267 else 2268 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 2269 & ' == Imposed occupation matrix' 2270 end if 2271 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 2272 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 2273 if (ndij==4) write(message,fmt='(4a)') trim(message)," for component ", & 2274 & trim(dspinc(ispden+2*(ndij/4)))," ==" 2275 call wrtout(std_out,message,wrt_mode) 2276 do ilm=1,2*lpawu+1 2277 if(cplex_dij==1)& 2278 & write(message,'(12(1x,9(1x,f10.5)))')& 2279 & (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1) 2280 if(cplex_dij==2)& 2281 & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 2282 & (paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden),jlm=1,2*lpawu+1) 2283 call wrtout(std_out,message,wrt_mode) 2284 end do 2285 end do 2286 2287 ! WARNING if symmetrization changes the matrix 2288 if (noccsym_error) then 2289 write(message, '(a,i4,6a)' ) & 2290 ' After symmetrization, imposed occupation matrix for atom ',iatom_tot,ch10,& 2291 & ' is different from dmatpawu value set in input file !',ch10,& 2292 & ' It is likely that dmatpawu does not match the symmetry operations of the system.',ch10,& 2293 & ' Action: change dmatpawu in input file or increase precision until 0.00001' 2294 ABI_WARNING(message) 2295 end if 2296 2297 end if ! impose_dmat/=0 2298 2299 ! ######################################################################################## 2300 ! # Rotate imposed occupation matrix in the non-diagonal basis 2301 ! ######################################################################################## 2302 if (usepawu/=0.and.impose_dmat/=0.and.dmatudiag_loc==2.and.(.not. cal_lmom)) then 2303 2304 lpawu=lcur;ldim=2*lpawu+1 2305 2306 ! Rotation of imposed nocc_mmp 2307 if (ndij/=4) then 2308 ABI_MALLOC(hdp2,(ldim,ldim)) 2309 do ispden=1,ndij 2310 call dgemm('n','t',ldim,ldim,ldim,one,& 2311 & paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim) 2312 call dgemm('n','n',ldim,ldim,ldim,one,& 2313 & noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim) 2314 end do ! ispden 2315 ABI_FREE(hdp2) 2316 else 2317 ABI_MALLOC(zhdp,(2*ldim,2*ldim)) 2318 ABI_MALLOC(zhdp2,(2*ldim,2*ldim)) 2319 do im2=1,ldim 2320 do im1=1,ldim 2321 zhdp( im1, im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1),zero,kind=dp) ! to be checked 2322 zhdp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2),zero,kind=dp) ! to be checked 2323 zhdp( im1,ldim+im2)=& 2324 & cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3),+paw_ij(iatom)%noccmmp(1,im1,im2,4),kind=dp) ! to be checked 2325 zhdp(ldim+im1, im2)=& 2326 & cmplx(paw_ij(iatom)%noccmmp(1,im2,im1,3),-paw_ij(iatom)%noccmmp(1,im2,im1,4),kind=dp) ! to be checked 2327 end do 2328 end do 2329 call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim) 2330 call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim) 2331 do jlm=1,ldim 2332 do ilm=1,ldim 2333 paw_ij(iatom)%noccmmp(1,ilm,jlm,1)= real(znoccmmp_tmp( ilm, jlm)) ! to be checked 2334 paw_ij(iatom)%noccmmp(1,ilm,jlm,2)= real(znoccmmp_tmp(ldim+ilm,ldim+jlm)) ! to be checked 2335 paw_ij(iatom)%noccmmp(1,ilm,jlm,3)= real(znoccmmp_tmp( ilm,ldim+jlm)) ! to be checked 2336 paw_ij(iatom)%noccmmp(1,ilm,jlm,4)=aimag(znoccmmp_tmp( ilm,ldim+jlm)) ! to be checked 2337 end do 2338 end do 2339 ABI_FREE(zhdp) 2340 ABI_FREE(zhdp2) 2341 end if 2342 2343 ! Printing of rotated imposed matrix 2344 do ispden=1,ndij 2345 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 2346 & ' == Imposed density matrix in original basis' 2347 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 2348 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 2349 if (ndij==4) write(message,fmt='(4a)') trim(message)," for component ", & 2350 & trim(dspin(ispden+2*(ndij/4)))," ==" 2351 call wrtout(std_out,message,wrt_mode) 2352 do ilm=1,2*lpawu+1 2353 write(message,'(12(1x,9(1x,f10.5)))') (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1) ! to be checked 2354 call wrtout(std_out,message,wrt_mode) 2355 end do 2356 end do ! ispden 2357 2358 end if ! dmatudiag_loc==2 2359 2360 if (usepawu/=0.and.dmatudiag_loc>0) then 2361 ABI_FREE(noccmmp_tmp) 2362 if (ndij==4) then 2363 ABI_FREE(znoccmmp_tmp) 2364 end if 2365 end if 2366 2367 paw_ij(iatom)%has_pawu_occ=2 2368 2369 end if ! lcur 2370 end do ! iatom 2371 2372 !Memory deallocation 2373 if (usepawu/=0.and.impose_dmat/=0) then 2374 do iatom_tot=1,natom 2375 lpawu=pawtab(typat(iatom_tot))%lpawu 2376 if (lpawu/=-1) then 2377 ABI_FREE(tmp_noccmmp(iatom_tot)%value) 2378 end if 2379 end do 2380 ABI_FREE(tmp_noccmmp) 2381 end if 2382 2383 !Destroy atom table used for parallelism 2384 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 2385 2386 DBG_EXIT("COLL") 2387 2388 end subroutine setnoccmmp
m_paw_correlations/setrhoijpbe0 [ Functions ]
[ Top ] [ m_paw_correlations ] [ Functions ]
NAME
setrhoijpbe0
FUNCTION
PAW local exact exchange only: Impose value of rhoij for f electrons using an auxiliairy file
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset initialized= if 0, the initialization of the gstate run is not yet finished istep=index of the number of steps in the routine scfcv istep_mix=index of the number of steps for the SCF mixing (can be <istep) mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms mpi_comm_read=MPI communicator containing all the processes reading the PBE0 file my_natom=number of atoms treated by current processor natom=number of atoms in cell ntypat=number of types of atoms in unit cell pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data typat(natom)=type integer for each atom in cell
SIDE EFFECTS
istep_mix=index of the number of steps for the SCF mixing (can be <istep) pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
NOTES
Only valid for f electrons !!!
SOURCE
2426 subroutine setrhoijpbe0(dtset,initialized,istep,istep_mix,& 2427 & mpi_comm_read,my_natom,natom,ntypat,pawrhoij,pawtab,typat,& 2428 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 2429 2430 !Arguments --------------------------------------------- 2431 !scalars 2432 integer,intent(in) :: initialized,istep,mpi_comm_read,my_natom,natom,ntypat 2433 integer,intent(inout) :: istep_mix 2434 integer,optional,intent(in) :: comm_atom 2435 type(dataset_type),intent(in) :: dtset 2436 !arrays 2437 integer,intent(in) :: typat(natom) 2438 integer,optional,target,intent(in) :: mpi_atmtab(:) 2439 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 2440 type(pawtab_type),intent(in) :: pawtab(ntypat) 2441 2442 !Local variables --------------------------------------- 2443 !scalars 2444 integer,parameter :: ll=3 2445 integer :: cplex_rhoij,iatom,iatom_tot,ierr,ii,ios,iread,irhoij,ispden,itypat,jj 2446 integer :: klmn,my_comm_atom,my_rank,nselect,nstep1,nstep1_abs,rhoijshft,rhoijsz 2447 logical :: my_atmtab_allocated,paral_atom,test0 2448 character(len=9),parameter :: filnam='rhoijpbe0' 2449 character(len=9),parameter :: dspin(6)=(/"up ","down ","up-up ","down-down","Re[up-dn]","Im[up-dn]"/) 2450 character(len=500) :: strg, message 2451 !arrays 2452 integer, allocatable :: nspden_tmp(:) 2453 integer,pointer :: my_atmtab(:) 2454 real(dp),allocatable :: rhoijtmp(:,:),rhoijtmp1(:,:),rhoijtmp2(:,:,:,:) 2455 2456 ! ********************************************************************* 2457 2458 DBG_ENTER("COLL") 2459 2460 !Some limitation 2461 if (my_natom>0) then 2462 if (pawrhoij(1)%qphase==2) then 2463 message='setrhoijpbe0 not compatible with qphase=2!' 2464 ABI_BUG(message) 2465 end if 2466 end if 2467 2468 !Test existence of file and open it 2469 inquire(file=filnam,iostat=ios,exist=test0) 2470 if(.not.test0) return 2471 2472 !Look for parallelisation over atomic sites 2473 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 2474 2475 !Test if exact-exch. is on f electrons 2476 test0=.false. 2477 do itypat=1,ntypat 2478 if (pawtab(itypat)%useexexch/=0.and.pawtab(itypat)%lexexch/=ll) test0=.true. 2479 end do 2480 if (test0) then 2481 write(message, '(3a,i1,a)' ) & 2482 & ' Local exact exchange: occ. matrix can only be imposed for l=',ll,' !' 2483 ABI_ERROR(message) 2484 end if 2485 2486 !============================================================ 2487 !===== First case: no parallelisation over atomic sites ===== 2488 !============================================================ 2489 2490 if (.not.paral_atom) then 2491 2492 ! Open file 2493 if (open_file(filnam,message,unit=77,form='formatted') /= 0) then 2494 ABI_ERROR(message) 2495 end if 2496 2497 ! Read step number and eventually exit 2498 nstep1=0;test0=.false. 2499 do while (.not.test0) 2500 read(77,'(A)') strg 2501 test0=(strg(1:1)/="#") 2502 if (test0) read(unit=strg,fmt=*) nstep1 2503 end do 2504 nstep1_abs=abs(nstep1) 2505 if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then 2506 close(77) 2507 ! Reinitalize mixing when rhoij is allowed to change; for experimental purpose... 2508 if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1 2509 return 2510 end if 2511 2512 ! Loop on atoms 2513 do iatom=1,natom 2514 itypat=typat(iatom) 2515 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 2516 2517 if (pawtab(itypat)%useexexch/=0) then 2518 2519 ! Set sizes depending on ll 2520 rhoijsz=4*ll+2 2521 rhoijshft=2*ll*ll 2522 2523 ! Uncompress rhoij 2524 ABI_MALLOC(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 2525 do ispden=1,pawrhoij(iatom)%nspden 2526 rhoijtmp=zero 2527 do irhoij=1,pawrhoij(iatom)%nrhoijsel 2528 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 2529 rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden) 2530 end do 2531 end do 2532 ! Read rhoij from file 2533 ABI_MALLOC(rhoijtmp1,(rhoijsz,rhoijsz)) 2534 do ispden=1,pawrhoij(iatom)%nspden 2535 do ii=1,rhoijsz 2536 test0=.false. 2537 do while (.not.test0) 2538 read(77,'(A)') strg 2539 test0=(strg(1:1)/="#") 2540 if (test0) read(unit=strg,fmt=*) (rhoijtmp1(ii,jj), jj=1,rhoijsz) 2541 end do 2542 end do 2543 2544 ! Impose rhoij 2545 do jj=1,rhoijsz 2546 do ii=1,jj 2547 rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp1(ii,jj) 2548 end do 2549 end do 2550 2551 end do 2552 ABI_FREE(rhoijtmp1) 2553 2554 ! Compress rhoij 2555 nselect=0 ; pawrhoij(iatom)%rhoijselect=0 2556 pawrhoij(iatom)%rhoijp=zero 2557 do klmn=1,pawrhoij(iatom)%lmn2_size 2558 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 2559 nselect=nselect+1 ; ii=cplex_rhoij*(nselect-1)+1 2560 do ispden=1,pawrhoij(iatom)%nspden 2561 pawrhoij(iatom)%rhoijp(ii,ispden)=rhoijtmp(klmn,ispden) 2562 end do 2563 pawrhoij(iatom)%rhoijselect(nselect)=klmn 2564 end if 2565 end do 2566 pawrhoij(iatom)%nrhoijsel=nselect 2567 ABI_FREE(rhoijtmp) 2568 2569 ! Print new rhoij 2570 do ispden=1,pawrhoij(iatom)%nspden 2571 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,& 2572 & ' == Imposed occupation matrix' 2573 if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 2574 if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 2575 if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)') trim(message)," for component ", & 2576 & trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," ==" 2577 call wrtout(std_out,message,'COLL') 2578 call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,& 2579 & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%lmn_size,ll,& 2580 & pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),& 2581 & 1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='COLL') 2582 end do 2583 2584 ! End loop on atoms 2585 end if 2586 end do 2587 2588 ! Close file 2589 close (77) 2590 2591 else 2592 2593 ! ============================================================ 2594 ! ====== 2nd case: no parallelisation over atomic sites ===== 2595 ! ============================================================ 2596 2597 my_rank=xmpi_comm_rank(mpi_comm_read) 2598 2599 ! Read step number and eventually exit 2600 iread=0 2601 if (my_rank==0) then 2602 if (open_file(filnam,message,unit=77,form='formatted') /=0 ) then 2603 ABI_ERROR(message) 2604 end if 2605 nstep1=0;test0=.false. 2606 do while (.not.test0) 2607 read(77,'(A)') strg 2608 test0=(strg(1:1)/="#") 2609 if (test0) read(unit=strg,fmt=*) nstep1 2610 end do 2611 nstep1_abs=abs(nstep1) 2612 if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then 2613 close(77) 2614 ! Reinitalize mixing when rhoij is allowed to change; for experimental purpose... 2615 if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1 2616 iread=1 2617 end if 2618 end if 2619 call xmpi_sum(iread,mpi_comm_read,ierr) 2620 if (iread/=0) return 2621 2622 ! Set up parallelism over atoms 2623 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 2624 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 2625 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 2626 2627 ! Store number of component for rhoij 2628 ABI_MALLOC(nspden_tmp,(natom)) 2629 nspden_tmp(:)=zero 2630 do iatom=1,my_natom 2631 iatom_tot=my_atmtab(iatom) 2632 nspden_tmp(iatom_tot)=pawrhoij(iatom)%nspden 2633 end do 2634 call xmpi_sum(nspden_tmp,mpi_comm_read,ierr) 2635 2636 ! To be improve if too much memory 2637 ABI_MALLOC(rhoijtmp2,(natom,4,rhoijsz,rhoijsz)) 2638 rhoijtmp2=zero 2639 2640 ! Read rhoij from file 2641 if (my_rank==0) then 2642 do iatom=1,natom 2643 itypat=typat(iatom) 2644 if (pawtab(itypat)%useexexch/=0) then 2645 rhoijsz=4*ll+2 2646 do ispden=1,nspden_tmp(iatom) 2647 do ii=1,rhoijsz 2648 test0=.false. 2649 do while (.not.test0) 2650 read(77,'(A)') strg 2651 test0=(strg(1:1)/="#") 2652 if (test0) read(unit=strg,fmt=*) (rhoijtmp2(iatom,ispden,ii,jj),jj=1,rhoijsz) 2653 end do 2654 end do 2655 end do 2656 end if 2657 end do 2658 end if 2659 call xmpi_sum(rhoijtmp2,mpi_comm_read,ierr) 2660 2661 ! Now, distribute rhoij 2662 do iatom=1,my_natom 2663 iatom_tot=my_atmtab(iatom) 2664 itypat=pawrhoij(iatom)%itypat 2665 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 2666 2667 if (pawtab(itypat)%useexexch/=0) then 2668 2669 ! Set sizes depending on ll 2670 rhoijsz=4*ll+2 2671 rhoijshft=2*ll*ll 2672 2673 ! Uncompress rhoij 2674 ABI_MALLOC(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 2675 do ispden=1,pawrhoij(iatom)%nspden 2676 rhoijtmp=zero 2677 do irhoij=1,pawrhoij(iatom)%nrhoijsel 2678 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 2679 rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden) 2680 end do 2681 2682 ! Impose rhoij 2683 do jj=1,rhoijsz 2684 do ii=1,jj 2685 rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp2(iatom_tot,ispden,ii,jj) 2686 end do 2687 end do 2688 2689 end do 2690 2691 ! Compress rhoij 2692 nselect=0 ; pawrhoij(iatom)%rhoijselect=0 2693 pawrhoij(iatom)%rhoijp=zero 2694 do klmn=1,pawrhoij(iatom)%lmn2_size 2695 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 2696 nselect=nselect+1 ; ii=cplex_rhoij*(nselect-1)+1 2697 do ispden=1,pawrhoij(iatom)%nspden 2698 pawrhoij(iatom)%rhoijp(ii,ispden)=rhoijtmp(klmn,ispden) 2699 end do 2700 pawrhoij(iatom)%rhoijselect(nselect)=klmn 2701 end if 2702 end do 2703 pawrhoij(iatom)%nrhoijsel=nselect 2704 ABI_FREE(rhoijtmp) 2705 2706 end if ! useexexch/=0 2707 2708 ! Print new rhoij 2709 do ispden=1,pawrhoij(iatom)%nspden 2710 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,' == Imposed occupation matrix' 2711 if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 2712 if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 2713 if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)') trim(message)," for component ", & 2714 & trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," ==" 2715 call wrtout(std_out,message,'PERS') 2716 call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,& 2717 & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%lmn_size,ll,& 2718 & pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),& 2719 & 1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='PERS') 2720 end do 2721 2722 ! end loop on atoms 2723 end do 2724 2725 ABI_FREE(nspden_tmp) 2726 ABI_FREE(rhoijtmp2) 2727 2728 ! Destroy atom table used for parallelism 2729 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 2730 2731 ! ============================================================ 2732 end if ! paral_atom 2733 2734 DBG_EXIT("COLL") 2735 2736 end subroutine setrhoijpbe0