Unstable behavior of the self-defined MPO's

DMRG and Numerical Methods
1

Dear Itensor Community,

I am currently working on structure factors of Bond Order Waves in 1D Hubbard model.
such a structure factor can be obtained through calculation of the following corelator:

\langle k_i k_j \rangle

,where k_i = \sum_{\sigma=(\uparrow,\downarrow)}c_{i+1,\sigma}c_{i,\sigma} + H.c .

Therefore one obtains a four point correlation function. Therefore, I have defined it in the following way:

function EHM_Correlators(psi::MPS,N::Int,sites)

sites = siteinds("Electron", N, conserve_qns=true);
replace_siteinds(psi,sites)

BCDW =  Array{Float64}(undef,N, N);
BSDW =  Array{Float64}(undef,N, N);
for i in 1:N-1
    for j in i+2:N-1

   # sites = siteinds("Electron", N, conserve_qns=true);
   # replace_siteinds(psi,sites)
        
    KCmpo = OpSum()
    KSmpo = OpSum()

    KCmpo += 1.0, "Cdagup", i+1,"Cup", i, "Cdagup", j+1,"Cup", j;
    KCmpo += 1.0, "Cdagup", i+1,"Cup", i, "Cdagdn", j+1,"Cdn", j;
    KCmpo += 1.0, "Cdagup", i+1,"Cup", i, "Cup", j+1,"Cdagup", j;
    KCmpo += 1.0, "Cdagup", i+1,"Cup", i, "Cdn", j+1,"Cdagdn", j;

    KCmpo += 1.0, "Cdagdn", i+1,"Cdn",i, "Cdagup", j+1,"Cup", j;
    KCmpo += 1.0, "Cdagdn", i+1,"Cdn",i , "Cdagdn", j+1,"Cdn", j;
    KCmpo += 1.0, "Cdagdn", i+1,"Cdn",i , "Cup", j+1,"Cdagup", j;
    KCmpo += 1.0, "Cdagdn", i+1,"Cdn",i , "Cdn", j+1,"Cdagdn", j;

    KCmpo += 1.0, "Cup", i+1,"Cdagup", i, "Cdagup", j+1,"Cup", j;
    KCmpo += 1.0, "Cup", i+1,"Cdagup", i, "Cdagdn", j+1,"Cdn", j;
    KCmpo += 1.0, "Cup", i+1,"Cdagup", i, "Cup", j+1,"Cdagup", j;
    KCmpo += 1.0, "Cup", i+1,"Cdagup", i, "Cdn", j+1,"Cdagdn", j;

    KCmpo += 1.0, "Cdn", i+1,"Cdagdn", i, "Cdagup", j+1,"Cup", j;
    KCmpo += 1.0, "Cdn", i+1,"Cdagdn", i, "Cdagdn", j+1,"Cdn", j;
    KCmpo += 1.0, "Cdn", i+1,"Cdagdn", i, "Cup", j+1,"Cdagup", j;
    KCmpo += 1.0, "Cdn", i+1,"Cdagdn", i, "Cdn", j+1,"Cdagdn", j;
    
    KC = MPO(KCmpo, sites);

    BCDW[i,j] = inner(psi',KC,psi)
    BCDW[j,i] = (inner(psi',KC,psi))'
    end
end

end

However, many times such matrix gives weird values, like 1e300, and other times it produces reasonable/correct results. I am kind of lost since cannot track the source of the problem further. I tried to replace site indeces inside the loop as well (as I commented in the code), since I thought the generation of MPO in every iteration might cause some problem with that.

I would very thankful for your help or any tips. If needed, I could also provide a minimal working example, but I quite sure the problem does not lie in the wave function.

2

I think where the error was. Let me confirm on monday.