Problems with CG_METHOD=0 [encountered in many different use cases]

beddalumia commented 2 years ago

What's up

As I was performing some initial tests on the incoming changes on fit_overhaul branch, I noticed something is wrong with the Weiss field analytical gradient, as it is currently on master. By wrong I mean that the fitted field is qualitatively different from the original one, so seriously wrong.

Note that this happens only if I request the analytic gradient: switching to CG_METHOD=1 (minimize routine) solves everything, switching only CG_GRAD=1 but remaining with the NR routine gives the correct fit but takes forever to compute, basically freezing for several minutes (✨holy minimize, mighty va10a✨).

Also note that the hybridization gradient appears to work fine, so maybe the actual problem lies within grad_g0and_replica. If that's the case I'll hit the wall very soon, when testing the gradients for the new 'elemental' norm. We'll see...

Evidence

<img src="https://user-images.githubusercontent.com/56808633/173280738-fdd20675-8c63-4a64-a8ab-6f3dd10bfc60.png" width="500" />

I did not include the real part in the figure for it is basically correct (but not totally: even with the numerical gradient I see some unwanted weight near the origin). Don't know if this could help, but I report it.
No restart file has been used, and here we are at one loop: we do not expect a brilliant fit, but a fair one yes.
I attach down here the input file, so that the problem could be reproduced with the current master head (dc046d84c8117d47c4e5fd623b43b03bf30316d7)

Notes

I've already added a warning for the user in commit b6351fff491015d0e0be7384fe04300803d8f493 (fit_overhaul), we might want to copy those lines on master without waiting for merge, if this is considered urgent.

Link to the lines

inputHM2D.conf

WMIXING=0.500000000                           !Mixing bath parameter
TS=2.500000000E-01                            !hopping parameter
NX=1                                          !Number of cluster sites in x direction
NY=1                                          !Number of cluster sites in y direction
NKX=30                                        !Number of kx point for BZ integration
NKY=30                                        !Number of ky point for BZ integration
NLAT=1                                        !Number of cluster sites
NORB=1                                        !Number of impurity orbitals (max 5).
NBATH=7                                       !Number of bath sites:(normal=>Nbath per orb)(hybrid=>Nbath total)(replica=>Nbath=Nreplica)
NSPIN=1                                       !Number of spin degeneracy (max 2)
ULOC=1.000000000,0.d0,0.d0,0.d0,0.d0          !Values of the local interaction per orbital (max 5)
UST=0.d0                                      !Value of the inter-orbital interaction term
JH=0.d0                                       !Hunds coupling
JX=0.d0                                       !S-E coupling
JP=0.d0                                       !P-H coupling
BETA=300.000000000                            !Inverse temperature, at T=0 is used as a IR cut-off.
XMU=0.d0                                      !Chemical potential. If HFMODE=T, xmu=0 indicates half-filling condition.
NLOOP=1                                       !Max number of DMFT iterations.
DMFT_ERROR=1.000000000E-05                    !Error threshold for DMFT convergence
SB_FIELD=1.000000000E-01                      !Value of a symmetry breaking field for magnetic solutions.
GF_FLAG=T                                     !flag to evaluate GFs and related quantities.
DM_FLAG=F                                     !flag to evaluate the cluster density matrix \rho_IMP = Tr_BATH(\rho))
ED_TWIN=F                                     !flag to reduce (T) or not (F,default) the number of visited sector using twin symmetry.
ED_SECTORS=F                                  !flag to reduce sector scan for the spectrum to specific sectors +/- ed_sectors_shift.
ED_SECTORS_SHIFT=1                            !shift to ed_sectors
ED_SPARSE_H=T                                 !flag to select  storage of sparse matrix H (mem--, cpu++) if TRUE, or direct on-the-fly H*v product (mem++, cpu--) if FALSE
ED_GF_SYMMETRIC=F                             !flag to assume Gij = Gji
ED_PRINT_SIGMA=T                              !flag to print impurity Self-energies
ED_PRINT_G=T                                  !flag to print impurity Greens function
ED_PRINT_G0=T                                 !flag to print non-interacting impurity Greens function
ED_VERBOSE=5                                  !Verbosity level: 0=almost nothing --> 5:all. Really: all
NSUCCESS=1                                    !Number of successive iterations below threshold for convergence
LMATS=5000                                    !Number of Matsubara frequencies.
LREAL=5000                                    !Number of real-axis frequencies.
LTAU=1024                                     !Number of imaginary time points.
LFIT=1000                                     !Number of Matsubara frequencies used in the \Chi2 fit.
NREAD=0.d0                                    !Objective density for fixed density calculations.
NERR=1.000000000E-04                          !Error threshold for fixed density calculations.
NDELTA=1.000000000E-01                        !Initial step for fixed density calculations.
NCOEFF=1.000000000                            !multiplier for the initial ndelta read from a file (ndelta-->ndelta*ncoeff).
WINI=-5.000000000                             !Smallest real-axis frequency
WFIN=5.000000000                              !Largest real-axis frequency
HFMODE=T                                      !Flag to set the Hartree form of the interaction (n-1/2). see xmu.
EPS=1.000000000E-02                           !Broadening on the real-axis.
CUTOFF=1.000000000E-09                        !Spectrum cut-off, used to determine the number states to be retained.
GS_THRESHOLD=1.000000000E-09                  !Energy threshold for ground state degeneracy loop up
HWBAND=2.000000000                            !half-bandwidth for the bath initialization: flat in -hwband:hwband
LANC_METHOD=arpack                            !select the lanczos method to be used in the determination of the spectrum. ARPACK (default), LANCZOS (T=0 only), DVDSON (no MPI)
LANC_NSTATES_SECTOR=2                         !Initial number of states per sector to be determined.
LANC_NSTATES_TOTAL=1                          !Initial number of total states to be determined.
LANC_NSTATES_STEP=2                           !Number of states added to the spectrum at each step.
LANC_NCV_FACTOR=10                            !Set the size of the block used in Lanczos-Arpack by multiplying the required Neigen (Ncv=lanc_ncv_factor*Neigen+lanc_ncv_add)
LANC_NCV_ADD=0                                !Adds up to the size of the block to prevent it to become too small (Ncv=lanc_ncv_factor*Neigen+lanc_ncv_add)
LANC_NITER=512                                !Number of Lanczos iteration in spectrum determination.
LANC_NGFITER=200                              !Number of Lanczos iteration in GF determination. Number of momenta.
LANC_TOLERANCE=1.000000000E-12                !Tolerance for the Lanczos iterations as used in Arpack and plain lanczos.
LANC_DIM_THRESHOLD=1024                       !Min dimension threshold to use Lanczos determination of the spectrum rather than Lapack based exact diagonalization.
CG_METHOD=0                                   !Conjugate-Gradient method: 0=NR, 1=minimize.
CG_GRAD=0                                     !Gradient evaluation method: 0=analytic (default), 1=numeric.
CG_FTOL=1.000000000E-05                       !Conjugate-Gradient tolerance.
CG_STOP=0                                     !Conjugate-Gradient stopping 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||).
CG_NITER=500                                  !Max. number of Conjugate-Gradient iterations.
CG_WEIGHT=1                                   !Conjugate-Gradient weight form: 1=1.0, 2=1/n , 3=1/w_n.
CG_SCHEME=weiss                               !Conjugate-Gradient fit scheme: delta or weiss.
CG_POW=2                                      !Fit power for the calculation of the Chi distance function as 1/L*|G0 - G0and|**cg_pow
CG_MINIMIZE_VER=F                             !Flag to pick old/.false. (Krauth) or new/.true. (Lichtenstein) version of the minimize CG routine
CG_MINIMIZE_HH=1.000000000E-04                !Unknown parameter used in the CG minimize procedure.
HFILE=hamiltonian                             !File where to retrieve/store the bath parameters.
IMPHFILE=inputHLOC.in                         !File read the input local H.
LOGFILE=6                                     !LOG unit.

beddalumia commented 2 years ago

NEWS (relative to fit_overhaul branch)

Commits 1bab32abf5295a43e58d54453268829bd249daf3 and 3240d4002ac08eef9bb3ee48e02b2779ebd88eba have shown that the issue is wider (hence the title change). Here I report some evidence and try to wrap a brief recap.

RECAP

With numerical gradients everything works, with both "Frobenius" and "Elemental" definition of $\chi^2$ distance.

	Deltaˆ	Weiss
Elemental	🟡	🟢
Frobenius	🟡	🟢

ˆThe Delta fits are tagged "yellow" for they have worse quality wrt the Weiss ones, if comparing with NR-CG results with DMFT_ED code (NR-CG freezes with CDMFT code, consistently. I don't know why). Yet minimize-CG results are all on the same level accross the two codes, and very similar to the "yellow" ones (so we can take them as "fairly good"). More info on the freezing in 3240d4002ac08eef9bb3ee48e02b2779ebd88eba commit message; plots of fitted functions are reported below. Note that instead minimize and NR give the very same results for the Weiss field, hence tagged green.

With analytical gradients we see some problems if we a) leave unchanged Frobenius gradients (wrt master branch) and b) port analytic gradients for elemental $\chi^2$ from current version in LIB_DMFT_ED. The situation is:

	Delta	Weiss
Elemental	🔴	🔴
Frobenius	🟡	🔴

If we swap sign in the elemental implementation of $\nabla\chi^2$ (thus diverging from LIB_DMFT_ED!) we get:

	Delta	Weiss
Elemental	🟡	🟡
Frobenius	🟡	🔴

If we further notice that commit 38bb300 had fixed the sign of \grad\chi^2(\Delta), but left untouched the sign of \grad\chi^2(g_0), so that applying the missing fix, we get:

	Delta	Weiss
Elemental	🟡	🟡
Frobenius	🟡	🟡

DETAILS

All plots with same input other than CG options, and at one loop. Solid Line: FG Dash-Dot: Fit

What I mean with "🔴" is

Delta	Weiss

What I mean with "🟡" is

Delta	Weiss

What I mean with "🟢" is

Deltaˆ	Weiss

ˆThis plot is the only one generated with DMFT_ED code, for this quality appears to be unreachable with minimize-CG, and NR-CG is de facto unavailable within CDMFT_ED.

beddalumia commented 2 years ago

Relevant update in SciFortran: https://github.com/QcmPlab/SciFortran/commit/c742471709f059fb66baede7bff5cf985470d9f9

Effects on this issue to be tested (it could probably fix the freezing with CG_METHOD=0 and CG_GRAD=1).

edit: it does not.

QcmPlab / CDMFT-LANC-ED