+1 vote
asked by (280 points)
edited by

I was wondering if you could help me with the following issue: For the sake of benchmarking my code, I am considering free fermions hopping on a 1D lattice, and set the filling to 2 fermions only. I furthermore set the spins of both electrons to be up. I use DMRG to solve this problem.
On the other hand, I solve it by diagonalizing the Hamiltonian matrix.

Then I look at correlation functions such as $\langle \psi | c^\daggeri cj | \psi \rangle$ in the above two cases, but the off-diagonal elements of such matrix differ in the two cases.

Could you please check and let me know how I can resolve this problem? I provided a minimal code below, sorry it got long. I used your own example of DMRG for the extended Hubbard model.

Thank you

using ITensors
using LinearAlgebra


N = 10
Npart = 2
t1 = 1.0

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

ampo = AutoMPO()
for b=1:N-1
ampo += -t1, "Cdagup", b,    "Cup", b+1
ampo += -t1, "Cdagup", b+1,  "Cup", b
ampo += -t1, "Cdagdn", b,    "Cdn", b+1
ampo += -t1, "Cdagdn", b+1,  "Cdn", b
end
H = MPO(ampo,sites)

sweeps = Sweeps(6)
maxdim!(sweeps,50,100,200,400,800,800)
cutoff!(sweeps,1E-12)

state = ["Emp" for n=1:N]
for i=1:Npart
    state[i] = "Up"
end

psi0 = randomMPS(sites, state, min(10, floor(Int,N/2) ))

# Check total number of particles:
@show flux(psi0)

# Start DMRG calculation:
energy, psi = dmrg(H, psi0, sweeps)



println("\nGround State Energy = $energy")


#Below, correlation functions using the DMRG solution are calculated:


cor_mat_mps = zeros(N,N);

for i = 1:N
    orthogonalize!(psi,i)

    s_ind_i = siteind(psi,i);

    for j =1:N
        if j < i
            continue
        end

        s_ind_j = siteind(psi,j);

        if i == j
            cor_mat_mps[i,i] = scalar( 
                psi[i] * prime( dag(psi[i]) , s_ind_i ) * op(s_ind_i,"Nup")
                );
        else

            prod_aux = psi[i];

            for n = i+1:j
                prod_aux *= psi[n];
            end

            cor_mat_mps[i,j]=   scalar(
                prod_aux * prime( dag(prod_aux) , s_ind_i, s_ind_j) *
                   op(s_ind_i,"Cdagup")*op(s_ind_j,"Cup")
                ); 

            cor_mat_mps[j,i]= scalar(
                prod_aux * prime( dag(prod_aux) , s_ind_i, s_ind_j) *
                   op(s_ind_j,"Cdagup")*op(s_ind_i,"Cup")
                ) ; 

        end
    end
end

# The same Hamiltonian is diagonalized below and correlation functions are calculated.

h_nonint = zeros(N,N);

for n = 1:N
    for m = 1:N
#         h_nonint[n,m] = 0.01*rand();
        if abs(m-n) == 1
            h_nonint[n,m] = -1. ;
        end

    end
end


diagonalization_energy = sum( eigen(Hermitian(h_nonint)).values[1:2] );

println("\nGround State Energy from diagonalization = $diagonalization_energy")


vecs = eigen(Hermitian(h_nonint)).vectors;

cor_direct = zeros( N,N );


for m = 1:N
    for n = 1:N
        cor_direct[m,n] = conj(vecs[m,1]) * vecs[n,1];
        cor_direct[m,n] += conj(vecs[m,2]) * vecs[n,2];
    end
end


println()

# below, the difference between correaltion functions calcualted in the two ways is shown, first only the diagonal terms and then all the correlation matrix.
println(  norm(diag(cor_direct-cor_mat_mps))  )
println(  norm(cor_direct-cor_mat_mps)  )

1 Answer

0 votes
answered by (70.1k points)
selected by
 
Best answer

Hi, so I see in your code that you are not including any Jordan-Wigner “F” string operators in your MPS correlation function calculations. This is necessary when computing pairs of operators such as C^dagger and C which have an odd fermion parity. (As a side comment, we are nearly done with a new feature that will remove this requirement, using a technique that makes tensor indices anticommutative, but it’s technical and still in the debugging phase.)

For some resources on how to do this and why, here is a sample code for the C++ version showing how to compute fermion correlators properly for the case of spinless fermions:
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/spinless_correlator_mps

For a discussion about Jordan-Wigner string and mappings to bosonic operators (all ITensors are bosonic in that sense), see this tutorial page: http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/fermions

Finally, the approach you are taking of making prod_aux where you combine many MPS tensors together looks like it will have an exponential cost, if I am understanding your code correctly there. This page has an explanation of the approach that is efficient for calculating MPS correlators: http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations

We can discuss some more if you have questions -

Best regards,
Miles

commented by (280 points)
Thank you for your elaborate response. I see the problem now; I was assuming that the Jordan-Wigner string is taken into account automatically when using C operators, since things worked properly with AutoMPO. But it seems it's not the case. Looking forward to working with this new feature that you mentioned.
Also, thanks for providing the more efficient approach to calculate the correlation functions, I was aware of this other approach but used the one above for clarity.
Thank you!
commented by (280 points)
I was wondering if you could help me with how I can adopt the code sample you provided above to the spinful case?
For the spinful case, when I have both spin up and down particles, the code does not give correct correlation function. Please see a sample code below.
Thanks

===============================
    using ITensors

    N = 10
    Npart = 2
    t1 = 1.0

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

    ampo = AutoMPO()
    for b=1:N-1
    ampo += -t1, "Cdagup", b,    "Cup", b+1
    ampo += -t1, "Cdagup", b+1,  "Cup", b
    ampo += -t1, "Cdagdn", b,    "Cdn", b+1
    ampo += -t1, "Cdagdn", b+1,  "Cdn", b
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(6)
    maxdim!(sweeps,50,100,200,400,800,800)
    cutoff!(sweeps,1E-12)

    state = ["Emp" for n=1:N]
    for i=1: floor(Int,Npart/2)
        state[i] = "Up"
        state[floor(Int,Npart/2)+i] = "Dn"
    end


    psi0 = randomMPS(sites, state, min(10, floor(Int,N/2) ))

    # Check total number of particles:
    @show flux(psi0)

    # Start DMRG calculation:
    energy, psi = dmrg(H, psi0, sweeps)



    println("\nGround State Energy = $energy")


    #Below, correlation functions using the DMRG solution are calculated:


    cor_mat_mps = zeros(N,N);

    for i = 1:N
        orthogonalize!(psi,i)

        s_ind_i = siteind(psi,i);
        l_ind_i = linkind(psi,i);

        for j =1:N
            if j < i
                continue
            end


            if i == j
                cor_mat_mps[i,i] = scalar(
                    psi[i] * prime( dag(psi[i]) , s_ind_i ) * op(s_ind_i,"Nup")
                    );
            else

                s_ind_j = siteind(psi,j);
                l_ind_j_left = linkind(psi,j-1);


                prod_aux = psi[i] * prime(dag(psi[i]), s_ind_i, l_ind_i);



                for n = i+1:j - 1
                    prod_aux *= psi[n] ;
                    prod_aux *= op(sites,"F",n) ;
                    prod_aux *= prime( dag( psi[n] ) ) ;
                end



                cor_mat_mps[i,j]=   scalar(
                    prod_aux * psi[j] * prime( dag( psi[j] ), s_ind_j, l_ind_j_left) *
                       op(s_ind_i,"Cdagup")*op(s_ind_j,"Cup")
                    );


                cor_mat_mps[j,i]= cor_mat_mps[i,j];

            end
        end
    end
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...