Numerical instabilities in Dyson Equation [KMH-MF]

beddalumia commented 2 years ago

Background

With commit cc0f13b961bc9d78ff984118515723b889ea6c35 we introduced an evaluation of "kinetic energy" $\langle H_0(k) \rangle$, exploiting the dmft_kinetic_energy(Hk,Sigma) function provided by DMFTtools. The strategy has been:

To build the noninteracting GFs by defining a "dummy" vanishing self-energy and feeding it, together with $H_0(k)$, to dmft_gloc_matsubara / dmft_gloc_realaxis.
To build the hartree-fock GFs by feeding again a vanishing self-energy to the same routines, but this time paired with $H_\mathrm{mf}(k)$ (since the interaction effects are treated at the mean-field level, they enter a corrected single-particle hamiltonian, and not the frequency dependent self-energy).
Compute the effective hartree-fock self-energy by means of a Dyson's equation, connecting the two Green's functions (i.e. formally acknowledge $H_\mathrm{mf}(k)$ to be interacting, hence transfer the mean-field correction from the single-particle hamiltonian, to a frequency independent self-energy).

Such a Dyson's equation reads:

$$\Sigma_\mathrm{mf}(z) = G0^{-1}(z) - G\mathrm{mf}^{-1}(z)$$

and has been implemented by:

function dyson_eq(G0,G) result(S)
    complex(8),dimension(Nlat,Nspin,Nspin,Norb,Norb,L) :: G0,G,S
    complex(8),dimension(Nlat,Nspin,Nspin,Norb,Norb)   :: Gnnn, G0nnn, Snnn
    complex(8),dimension(Nlso,Nlso)                    :: Glso, G0lso, Slso
    !
    do i = 1,L
       G0nnn = G0(:,:,:,:,:,i) ; Gnnn = G(:,:,:,:,:,i)
       G0lso = nnn2lso(G0nnn)  ; Glso = nnn2lso(Gnnn)
       call inv(G0lso)         ; call inv(Glso)       
       Slso = G0lso - Glso     ; Snnn = lso2nnn(Slso)
       S(:,:,:,:,:,i) = Snnn
    enddo
    !
 end function dyson_eq

The resulting kinetic energy appeared to be "good enough", if compared to the analogous DMFT result:

so that we called it a day.

The problem

But inspecting the actual printed self-energies we recognize something is affecting seriously our computation: we expect frequency independent self-energies, we get a total mess.

And by total mess I mean:

Real Axis	Matsubara

Arguably the Matsubara axis is quite a lot less crazy, but still:

it's not even causal
for sure not really constant

First Guess

I tried to understand if we have a problem on how we implement the Dyson's equation by re-implementing it in MATLAB, by reading from file all the components of $G0(z), G\mathrm{mf}(z)$, by means of the following script:

import plotDMFT.*

G011 = spectral_load('G0mats_l11_s1_iw__indx000001.dat');  
G012 = spectral_load('G0mats_l11_s1_iw__indx000002.dat');
G021 = spectral_load('G0mats_l11_s2_iw__indx000001.dat');
G022 = spectral_load('G0mats_l11_s2_iw__indx000002.dat');

G11 = spectral_load('Gmats_l11_s1_iw__indx000001.dat');
G12 = spectral_load('Gmats_l11_s1_iw__indx000002.dat');
G21 = spectral_load('Gmats_l11_s2_iw__indx000001.dat');
G22 = spectral_load('Gmats_l11_s2_iw__indx000002.dat');

S11 = spectral_load('Smats_l11_s1_iw__indx000001.dat');
S12 = spectral_load('Smats_l11_s1_iw__indx000002.dat');
S21 = spectral_load('Smats_l11_s2_iw__indx000001.dat');
S22 = spectral_load('Smats_l11_s2_iw__indx000002.dat');

domain = G11.zeta;

G0_matrix = zeros(2);
G_matrix  = zeros(2);
S_matrix  = zeros(2);

for iw = 1:length(domain)

    fprintf('iw = %f\n\n',domain(iw));

    G0_matrix(1,1) = G011.real(iw) + 1j * G011.imag(iw);
    G0_matrix(1,2) = G012.real(iw) + 1j * G012.imag(iw);
    G0_matrix(2,1) = G021.real(iw) + 1j * G021.imag(iw);
    G0_matrix(2,2) = G022.real(iw) + 1j * G022.imag(iw);

    G_matrix(1,1) = G11.real(iw) + 1j * G11.imag(iw);
    G_matrix(1,2) = G12.real(iw) + 1j * G12.imag(iw);
    G_matrix(2,1) = G21.real(iw) + 1j * G21.imag(iw);
    G_matrix(2,2) = G22.real(iw) + 1j * G22.imag(iw);

    S_matrix(1,1) = S11.real(iw) + 1j * S11.imag(iw);
    S_matrix(1,2) = S12.real(iw) + 1j * S12.imag(iw);
    S_matrix(2,1) = S21.real(iw) + 1j * S21.imag(iw);
    S_matrix(2,2) = S22.real(iw) + 1j * S22.imag(iw);

    D_matrix = inv(G0_matrix) - inv(G_matrix);

    if not(isequal(D_matrix,S_matrix))
        disp(D_matrix-S_matrix)
    else
        disp('ok!')
    end

end  % tried with real-axis too: same warning, similar pop-up frequency

And actually got instantly a quite unsettling warning:

>> Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = xxxxxxxx.

where typical values of RCOND fell in the (E-20,E-16) range.

SciFortran's inv() does not complain at all, but I do not see any reason for the situation to be different across the two languages. Most probably the problems lies in the zeros, for real frequencies, and the exponentially vanishing tail, for the Matsubara axis.

Shame on me to not foresee such an obvious numerical pitfall.

beddalumia commented 2 years ago

Notes:

Probably the kinetic energy is not so much wrong just because the real part of the Matsubara self-energy is almost constant and maybe near the right value.
For the in-plane magnetization actually we expect all the spin-diagonal components of the self-energy to vanish, and currently we do not print off-diagonal ones. We should change the iprint value from 4 to [?].
The fixing plan is to build the self-energy explicitly, i.e. from the order parameters. The actual formula? I have to think...

beddalumia commented 2 years ago

Idea for the explicit construction of ∑(ω)

We know that:

$\Sigma(\omega)$ has to be a constant at the Hartree-Fock level. Hence $\Sigma\mathrm{mf}(\infty) = \Sigma\mathrm{mf}(0) = \Sigma_\mathrm{mf}(\omega),\enspace\forall\omega$.
$\Sigma(0)$ can be used to renormalize $H0(k)$, preserving the topology of the interacting model, as $H\mathrm{top}(k) = H_0(k) + \Sigma(0)$
The Hartree-Fock hamiltonian is just a static renormalization of the noninteracting one, as $H_\mathrm{mf}(k) = H0(k) + \delta H\mathrm{mf}$.
$H_\mathrm{mf}(k)$ is de facto a single-particle hamiltonian, so that its topological hamiltonian should just be itself.

All of this should imply that:

$$H\mathrm{top} = H\mathrm{mf} \quad\Longrightarrow\quad \Sigma\mathrm{mf} = \delta H\mathrm{mf}$$

All this is probably too much of a convoluted reasoning for the matter, but should be correct afterall.

beddalumia commented 2 years ago

Can we really extend the reasoning to the imaginary axis?

i.e. can we state that $\Sigma(z)$ is constant on the whole complex plane?

beddalumia commented 2 years ago

Apparently yes... we can assume $\Sigma(z) = \delta H_\mathrm{mf}$ on the whole complex plane and we get excellent results.

Kinetic energy

We get same qualitative behavior, but improved comparison with the DMFT computation:

How could it be so good before?

Almost unbelievably we had, on the imaginary axis, a noncausal, nonconstant self-energy (so mightily wrong), but actually slow-varying only near the origin, and most impressively, almost correct at mid to large frequencies. So, in a cursed numerical sense, we were somewhat near the right solution, as far as the kinetic energy computation was concerned.

:sparkles:Marvelous!:sparkles:

beddalumia commented 2 years ago

IMPORTANT NOTE TO SELF

Most probably the true reason for the miserable failure of the Dyson approach lies in the actual non validity of such an equation relating $G0^\mathrm{loc}, \enspace G\mathrm{mf}^\mathrm{loc}, \enspace\mathrm{and} \enspace \Sigma_\mathrm{mf}^\mathrm{loc}$.

In fact the Dyson's equation that for sure olds is:

$$\Sigma_\mathrm{mf}(k,z) = G0^{-1}(k,z) - G\mathrm{mf}^{-1}(k,z)$$

which is indeed equiivalent to:

$$\sumk \Sigma\mathrm{mf}(k,z) = \sum_k G_0^{-1}(k,z) - \sumk G\mathrm{mf}^{-1}(k,z)$$

but not to:

$$\Sigma_\mathrm{mf}^\mathrm{loc}(k,z) = (G0^\mathrm{loc})^{-1}(k,z) - (G\mathrm{mf}^\mathrm{loc})^{-1}(k,z)$$

for that we cannot in any way move the sum over momentum into the matrix (or even the scalar) inversion operator.

beddalumia / KMHproject