Skip to content

Commit

Permalink
Merge branch 'master' of github.com:aamaricci/EDIpack2.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@aamaricci
aamaricci committed May 21, 2024
2 parents 26154a7 + 7d80916 commit c5235d4
Show file tree
Hide file tree
Showing 8 changed files with 320 additions and 92 deletions.
2 changes: 1 addition & 1 deletion src/ED_BATH/ED_BATH_AUX.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,7 @@ function Hgeneral_build(lambdavec) result(H)
complex(8),dimension(Nnambu*Nspin,Nnambu*Nspin,Norb,Norb) :: H
!
if(.not.Hgeneral_status)STOP "ERROR Hgeneral_build: Hgeneral_basis is not setup"
allocate(lambda(size(Hreplica_basis)));lambda=1d0
allocate(lambda(size(Hgeneral_basis)));lambda=1d0
if(present(lambdavec))then
if(size(lambdavec)/=size(Hgeneral_basis)) STOP "ERROR Hgeneral_build: Wrong coefficient vector size"
lambda = lambdavec
Expand Down
2 changes: 1 addition & 1 deletion src/ED_BATH/ED_BATH_DMFT.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -462,7 +462,7 @@ subroutine write_dmft_bath(dmft_bath_,unit)
"Vbk_l"//reg(str(iorb))//"_s"//reg(str(ispin)),&
iorb=1,Norb), ispin=1,Nspin)
do i=1,Nbath
write(unit,"(90(ES21.12,1X))")((&
write(unit_,"(90(ES21.12,1X))")((&
dmft_bath_%e(ispin,iorb,i),&
dmft_bath_%v(ispin,iorb,i),&
dmft_bath_%u(ispin,iorb,i),&
Expand Down
254 changes: 240 additions & 14 deletions src/ED_BATH/ED_FIT_GENERAL.f90
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
MODULE ED_FIT_GENERAL
USE ED_FIT_COMMON
USE SF_SPIN, only: pauli_sigma_z
USE SF_LINALG, only: kron, diag
USE SF_LINALG, only: kron, diag, trace

implicit none
private
Expand All @@ -10,7 +10,9 @@ MODULE ED_FIT_GENERAL
public :: chi2_fitgf_general
public :: chi2_fitgf_general_superc


complex(8),dimension(:,:,:,:,:),allocatable :: FGmatrix


contains


Expand Down Expand Up @@ -43,6 +45,12 @@ subroutine chi2_fitgf_general(fg,bath_)
check= check_bath_dimension(bath_)
if(.not.check)stop "chi2_fitgf_general error: wrong bath dimensions"
!
if(cg_pow/=2.AND.cg_norm=="frobenius")then
print *, "WARNING: CG_POW must be 2 for a meaningful definition of the Frobenius norm."
print *, " we'll still let you go ahead with the desired input, but please be "
print *, " be aware that CG_POW is not doing what you would expect for a chi^q"
endif
!
call allocate_dmft_bath(dmft_bath)
call set_dmft_bath(bath_,dmft_bath)
allocate(array_bath(size(bath_)-1))
Expand Down Expand Up @@ -78,6 +86,7 @@ subroutine chi2_fitgf_general(fg,bath_)
!
Ldelta = Lfit ; if(Ldelta>size(fg,5))Ldelta=size(fg,5)
!
allocate(FGmatrix(Nspin,Nspin,Norb,Norb,Ldelta))
allocate(Gdelta(totNso,Ldelta))
allocate(Xdelta(Ldelta))
allocate(Wdelta(Ldelta))
Expand All @@ -102,8 +111,11 @@ subroutine chi2_fitgf_general(fg,bath_)
if(ed_verbose>3)write(Logfile,"(A)")&
"DEBUG chi2_fitgf_general: cg_method:"//str(cg_method)//&
", cg_grad:"//str(cg_grad)//&
", cg_scheme:"//str(cg_scheme)
", cg_scheme:"//str(cg_scheme)//&
", cg_norm:"//str(cg_norm)
#endif
!
FGmatrix=fg
!
select case(cg_method) !0=NR-CG[default]; 1=CG-MINIMIZE
case default
Expand Down Expand Up @@ -249,7 +261,23 @@ end subroutine chi2_fitgf_general
!+-----------------------------------------------------------------+
!PURPOSE: Evaluate the \chi^2 distance of \Delta_Anderson function.
!+-----------------------------------------------------------------+
function chi2_delta_general(a) result(chi2)
function chi2_delta_general(a) result(chi2)
real(8),dimension(:) :: a
real(8) :: chi2
!
select case(cg_norm)
case ("elemental")
chi2 = chi2_delta_general_elemental(a)
case ("frobenius")
chi2 = chi2_delta_general_frobenius(a)
case default
stop "chi2_fitgf_general error: cg_norm != [elemental,frobenius]"
end select
!
end function chi2_delta_general


function chi2_delta_general_elemental(a) result(chi2)
real(8),dimension(:) :: a
real(8) :: chi2
real(8),dimension(totNso) :: chi2_so
Expand Down Expand Up @@ -279,13 +307,57 @@ function chi2_delta_general(a) result(chi2)
if(ed_verbose>3)write(Logfile,"(A,ES10.2)")"DEBUG chi2_delta_general. Chi**2:",chi2
#endif
!
end function chi2_delta_general
end function chi2_delta_general_elemental
!

!> FROBENIUS NORM: global \chi^2 for all components, only i\omega are weighted
function chi2_delta_general_frobenius(a) result(chi2)
real(8),dimension(:) :: a
real(8) :: chi2
real(8),dimension(Ldelta) :: chi2_freq
complex(8),dimension(Nspin,Nspin,Norb,Norb,Ldelta) :: Delta
complex(8),dimension(Nspin*Norb,Nspin*Norb) :: Delta_so
integer :: l
!
#ifdef _DEBUG
if(ed_verbose>5)write(Logfile,"(A,"//str(size(a))//"ES10.2)")"DEBUG chi2_delta_general_frobenius. a:",a
#endif
!
Delta = delta_general(a)
!
do l=1,Ldelta
Delta_so = nn2so_reshape(delta(:,:,:,:,l) - FGmatrix(:,:,:,:,l),Nspin,Norb)
chi2_freq(l) = sqrt(trace(matmul(Delta_so,conjg(transpose(Delta_so)))))
enddo
!
chi2 = sum(chi2_freq**cg_pow/Wdelta) !Weighted sum over matsubara frqs
chi2 = chi2/Ldelta/(Nspin*Norb) !Normalization over {iw} and Nlso
!
#ifdef _DEBUG
if(ed_verbose>3)write(Logfile,"(A,ES10.2)")"DEBUG chi2_delta_general_frobenius. chi2:",chi2
#endif
!
end function chi2_delta_general_frobenius
!
!+--------------------------------------------------------------+
!PURPOSE: Evaluate the \chi^2 distance of G_0_Anderson function.
!+--------------------------------------------------------------+
function chi2_weiss_general(a) result(chi2)
function chi2_weiss_general(a) result(chi2)
real(8),dimension(:) :: a
real(8) :: chi2
!
select case(cg_norm)
case ("elemental")
chi2 = chi2_weiss_general_elemental(a)
case ("frobenius")
chi2 = chi2_weiss_general_frobenius(a)
case default
stop "chi2_fitgf_general error: cg_norm != [elemental,frobenius]"
end select
!
end function chi2_weiss_general

function chi2_weiss_general_elemental(a) result(chi2)
real(8),dimension(:) :: a
real(8) :: chi2
real(8),dimension(totNso) :: chi2_so
Expand Down Expand Up @@ -315,21 +387,61 @@ function chi2_weiss_general(a) result(chi2)
if(ed_verbose>3)write(Logfile,"(A,ES10.2)")"DEBUG chi2_weiss_general. Chi**2:",chi2
#endif
!
end function chi2_weiss_general

end function chi2_weiss_general_elemental

!> FROBENIUS NORM: global \chi^2 for all components, only i\omega are weighted
function chi2_weiss_general_frobenius(a) result(chi2)
real(8),dimension(:) :: a
real(8) :: chi2
real(8),dimension(Ldelta) :: chi2_freq
complex(8),dimension(Nspin,Nspin,Norb,Norb,Ldelta) :: g0and
complex(8),dimension(Nspin*Norb,Nspin*Norb) :: Delta_so
integer :: l
!
#ifdef _DEBUG
if(ed_verbose>5)write(Logfile,"(A,"//str(size(a))//"ES10.2)")"DEBUG chi2_weiss_general_frobenius. a:",a
#endif
!
g0and = g0and_general(a)
!
do l=1,Ldelta
Delta_so = nn2so_reshape(g0and(:,:,:,:,l) - FGmatrix(:,:,:,:,l),Nspin,Norb)
chi2_freq(l) = sqrt(trace(matmul(Delta_so,conjg(transpose(Delta_so)))))
enddo
!
chi2 = sum(chi2_freq**cg_pow/Wdelta) !Weighted sum over matsubara frqs
chi2 = chi2/Ldelta/(Nspin*Norb) !Normalization over {iw} and Nlso
!
#ifdef _DEBUG
if(ed_verbose>3)write(Logfile,"(A,ES10.2)")"DEBUG chi2_weiss_general_frobenius. chi2:",chi2
#endif
!
end function chi2_weiss_general_frobenius



!######################################################################
! THESE PROCEDURES EVALUATE THE >GRADIENTS< OF THE \chi^2 TO MINIMIZE.
!######################################################################
!
function grad_chi2_delta_general(a) result(dchi2)
real(8),dimension(:) :: a
real(8),dimension(size(a)) :: dchi2
!
select case(cg_norm)
case ("elemental")
dchi2 = grad_chi2_delta_general_elemental(a)
case ("frobenius")
dchi2 = grad_chi2_delta_general_frobenius(a)
case default
stop "chi2_fitgf_general error: cg_norm != [elemental,frobenius]"
end select
!
end function grad_chi2_delta_general
!
!+---------------------------------------------------------------------+
!PURPOSE: Evaluate the gradient \Grad\chi^2 of \Delta_Anderson function.
!+---------------------------------------------------------------------+
function grad_chi2_delta_general(a) result(dchi2)
function grad_chi2_delta_general_elemental(a) result(dchi2)
real(8),dimension(:) :: a
real(8),dimension(size(a)) :: dchi2
real(8),dimension(totNso,size(a)) :: df
Expand Down Expand Up @@ -366,14 +478,73 @@ function grad_chi2_delta_general(a) result(dchi2)
if(ed_verbose>4)write(Logfile,"(A,"//str(size(a))//"ES10.2)")"DEBUG grad_chi2_delta_general. dChi**2:",dchi2
#endif
!
end function grad_chi2_delta_general
!
end function grad_chi2_delta_general_elemental
!
!> FROBENIUS NORM: global \chi^2 for all components, only i\omega are weighted
function grad_chi2_delta_general_frobenius(a) result(dchi2)
real(8),dimension(:) :: a
real(8),dimension(size(a)) :: dchi2
real(8),dimension(Ldelta,size(a)) :: df
complex(8),dimension(Nspin,Nspin,Norb,Norb,Ldelta) :: Delta
complex(8),dimension(Nspin,Nspin,Norb,Norb,Ldelta,size(a)) :: dDelta
complex(8),dimension(Ldelta) :: Ftmp
real(8),dimension(Ldelta,size(a)) :: dChi_freq
integer :: i,j,idelta,iorb,jorb,ispin,jspin
!
Delta = delta_general(a)
dDelta = grad_delta_general(a)
Ftmp=zero
df=zero
!
do idelta=1,Ldelta
do ispin=1,Nspin
do jspin=1,Nspin
do iorb=1,Norb
do jorb=1,Norb
!
Ftmp(idelta) = Ftmp(idelta) + abs(Delta(ispin,jspin,iorb,jorb,idelta)-FGmatrix(ispin,jspin,iorb,jorb,idelta))**2
do j=1,size(a)
df(idelta,j) = df(idelta,j) + &
real(Delta(ispin,jspin,iorb,jorb,idelta) - FGmatrix(ispin,jspin,iorb,jorb,idelta)) * &
real(dDelta(ispin,jspin,iorb,jorb,idelta,j)) + &
imag(Delta(ispin,jspin,iorb,jorb,idelta) - FGmatrix(ispin,jspin,iorb,jorb,idelta)) * &
imag(dDelta(ispin,jspin,iorb,jorb,idelta,j))
enddo
enddo
enddo
enddo
enddo
Ftmp(idelta) = cg_pow * (sqrt(Ftmp(idelta))**(cg_pow-2)) / Wdelta(idelta)
dchi_freq(idelta,:) = Ftmp(idelta) * df(idelta,:)
enddo
!
dchi2 = sum(dchi_freq,1)/Ldelta/(Nspin*Norb)
!
#ifdef _DEBUG
if(ed_verbose>4)write(Logfile,"(A,"//str(size(a))//"ES10.2)")"DEBUG grad_chi2_delta_general_frobenius. dChi2:",dchi2
#endif
!
end function grad_chi2_delta_general_frobenius
!
!+------------------------------------------------------------------+
!PURPOSE: Evaluate the gradient \Grad\chi^2 of G0_Anderson function.
!+------------------------------------------------------------------+
function grad_chi2_weiss_general(a) result(dchi2)
function grad_chi2_weiss_general(a) result(dchi2)
real(8),dimension(:) :: a
real(8),dimension(size(a)) :: dchi2
!
select case(cg_norm)
case ("elemental")
dchi2 = grad_chi2_weiss_general_elemental(a)
case ("frobenius")
dchi2 = grad_chi2_weiss_general_frobenius(a)
case default
stop "chi2_fitgf_general error: cg_norm != [elemental,frobenius]"
end select
!
end function grad_chi2_weiss_general

function grad_chi2_weiss_general_elemental(a) result(dchi2)
real(8),dimension(:) :: a
real(8),dimension(size(a)) :: dchi2
real(8),dimension(totNso,size(a)) :: df
Expand Down Expand Up @@ -410,8 +581,63 @@ function grad_chi2_weiss_general(a) result(dchi2)
if(ed_verbose>4)write(Logfile,"(A,"//str(size(a))//"ES10.2)")"DEBUG grad_chi2_weiss_general. dChi**2:",dchi2
#endif
!
end function grad_chi2_weiss_general
end function grad_chi2_weiss_general_elemental

!> FROBENIUS NORM: global \chi^2 for all components, only i\omega are weighted
function grad_chi2_weiss_general_frobenius(a) result(dchi2)
real(8),dimension(:) :: a
real(8),dimension(size(a)) :: dchi2
real(8),dimension(Ldelta,size(a)) :: df
complex(8),dimension(Nspin,Nspin,Norb,Norb,Ldelta) :: g0and
complex(8),dimension(Nspin,Nspin,Norb,Norb,Ldelta,size(a)) :: dg0and
complex(8),dimension(Ldelta) :: Ftmp
real(8),dimension(Ldelta,size(a)) :: dChi_freq
integer :: i,j,idelta,iorb,jorb,ispin,jspin
!
!
print*, " "
print*, "WARNING: the analytic gradient of the Weiss field Frobenius distance "
print*, " has been found giving /QUALITATIVELY WRONG/ fitted spectra! "
print*, " > the issue has emerged in some Nlat=Nspin=Norb=1 test runs "
print*, " > while the bug is investigated please switch cg_scheme to "
print*, " 'delta' or cg_norm to 'elemental' (or give up analytic cg)"
print*, " "
!
!
g0and = g0and_general(a)
dg0and = grad_g0and_general(a)
Ftmp=zero
df=zero
!
do idelta=1,Ldelta
do ispin=1,Nspin
do jspin=1,Nspin
do iorb=1,Norb
do jorb=1,Norb
!
Ftmp(idelta) = Ftmp(idelta) + abs(g0and(ispin,jspin,iorb,jorb,idelta)-FGmatrix(ispin,jspin,iorb,jorb,idelta))**2
do j=1,size(a)
df(idelta,j) = df(idelta,j) + &
real(g0and(ispin,jspin,iorb,jorb,idelta) - FGmatrix(ispin,jspin,iorb,jorb,idelta)) * &
real(dg0and(ispin,jspin,iorb,jorb,idelta,j)) + &
imag(g0and(ispin,jspin,iorb,jorb,idelta) - FGmatrix(ispin,jspin,iorb,jorb,idelta)) * &
imag(dg0and(ispin,jspin,iorb,jorb,idelta,j))
enddo
enddo
enddo
enddo
enddo
Ftmp(idelta) = cg_pow * (sqrt(Ftmp(idelta))**(cg_pow-2)) / Wdelta(idelta)
dchi_freq(idelta,:) = Ftmp(idelta) * df(idelta,:)
enddo
!
dchi2 = sum(dchi_freq,1)/Ldelta/(Nspin*Norb)
!
#ifdef _DEBUG
if(ed_verbose>4)write(Logfile,"(A,"//str(size(a))//"ES10.2)")"DEBUG grad_chi2_weiss_general_frobenius. dChi2:",dchi2
#endif
!
end function grad_chi2_weiss_general_frobenius



Expand Down
2 changes: 2 additions & 0 deletions src/ED_INPUT_VARS.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@ MODULE ED_INPUT_VARS
integer :: cg_stop !fit stop condition:0-3, 0=C1.AND.C2, 1=C1, 2=C2 with C1=|F_n-1 -F_n|<tol*(1+F_n), C2=||x_n-1 -x_n||<tol*(1+||x_n||).
integer :: cg_Weight !CGfit mode 0=1, 1=1/n , 2=1/w_n weight
integer :: cg_pow !fit power to generalize the distance as |G0 - G0and|**cg_pow
character(len=12) :: cg_norm !frobenius/elemental (for now only in general bath)
logical :: cg_minimize_ver !flag to pick old (Krauth) or new (Lichtenstein) version of the minimize CG routine
real(8) :: cg_minimize_hh !unknown parameter used in the CG minimize procedure.
!
Expand Down Expand Up @@ -248,6 +249,7 @@ subroutine ed_read_input(INPUTunit)
call parse_input_variable(cg_niter,"CG_NITER",INPUTunit,default=500,comment="Max. number of Conjugate-Gradient iterations.")
call parse_input_variable(cg_weight,"CG_WEIGHT",INPUTunit,default=1,comment="Conjugate-Gradient weight form: 1=1.0, 2=1/n , 3=1/w_n.")
call parse_input_variable(cg_scheme,"CG_SCHEME",INPUTunit,default='weiss',comment="Conjugate-Gradient fit scheme: delta or weiss.")
call parse_input_variable(cg_norm,"CG_NORM",INPUTunit,default='elemental',comment="Conjugate-Gradient norm definition: elemental (default) or frobenius.")
call parse_input_variable(cg_pow,"CG_POW",INPUTunit,default=2,comment="Fit power for the calculation of the generalized distance as |G0 - G0and|**cg_pow")
call parse_input_variable(cg_minimize_ver,"CG_MINIMIZE_VER",INPUTunit,default=.false.,comment="Flag to pick old/.false. (Krauth) or new/.true. (Lichtenstein) version of the minimize CG routine")
call parse_input_variable(cg_minimize_hh,"CG_MINIMIZE_HH",INPUTunit,default=1d-4,comment="Unknown parameter used in the CG minimize procedure.")
Expand Down
Loading

0 comments on commit c5235d4

Please sign in to comment.