Sz(M)Sz(i) Correlation Function in a Thermal State

Sambu_1234 · December 2, 2024, 6:50pm

Hi everyone,
I want to calculate the S_{z}(M)S_{z}(i) correlation function in the thermal state for a Heisenberg spin-1/2 model. To construct the thermal state, I used the ancilla method, where each pair of sites (physical + ancilla) was initialized in the state 1/ sqrt 2 (∣up up⟩+∣down down⟩), then performed imaginary time evolution up to \beta/2 using TDVP method. Odd sites are the physical sites, and even sites are the ancilla sites.

Here is the code attached:

using namespace itensor;

int main()
{
int N = 10;
int num_sites = 2 * N;
auto beta = 0.1;
auto nsteps = 5;
auto delta_beta = 0.05; // Increment for beta
auto num_beta_steps = 5;

auto sites = SpinHalf(num_sites, {"ConserveQNs=", false}); // 2N sites

auto psi = MPS(sites);
println(psi);
for(int n = 1; n <= 2*N;  n += 2)
{
    auto s1 = sites(n);
    auto s2 = sites(n+1);
    auto wf = ITensor(s1,s2); 
    wf.set(s1(1),s2(1), 1.0 / std::sqrt(2.0));
    wf.set(s1(2),s2(2), 1.0 / std::sqrt(2.0));
    ITensor D; 
    psi.ref(n) = ITensor(s1);
    psi.ref(n+1) = ITensor(s2);
    svd(wf,psi.ref(n),D,psi.ref(n+1));
    psi.ref(n) *= D;
}
Print(psi);

psi.position(1);

auto ampo = AutoMPO(sites);

for (int j = 1; j < num_sites-1; j += 2) {
    ampo += 0.5, "S+", j, "S-", j + 2;  
    ampo += 0.5, "S-", j, "S+", j + 2; 
    ampo += 1.0, "Sz", j, "Sz", j + 2; 
}

auto H = toMPO(ampo);

// TDVP parameters
auto sweeps = Sweeps(1);
sweeps.maxdim() = 64;
sweeps.cutoff() = 1E-10;
sweeps.niter() = 10;

println("----------------------------------------GSE-TDVP---------------------------------------");

auto current_beta = beta;
for (int bstep = 0; bstep < num_beta_steps; ++bstep)
{
    auto psi1 = psi; // Initialize or reuse the MPS state
    current_beta = beta + bstep * delta_beta;
    for (int n = 1; n <= nsteps; ++n)
    {
     psi1 = tdvp(psi1, H, -(current_beta) / (2 * nsteps), sweeps, {
            "truncate", true,
            "DoNormalize", true,
            "Quiet", true,
            "NumCenter", 2});

    Print(psi1);

// m = middle spin of the chain (only physical sites considered)

    int m = 0.5*(num_sites-2);

   auto op_m = op(sites, "Sz", m);

for (auto j = 1; j <= (num_sites-1); j += 2)
{
    // Make the operator for site j
    auto op_j = op(sites, "Sz", j);

        println("Correlation (<S_z(M) S_z_" ">) = ", m," ", j);

    if (m < j)
    {
        psi1.position(m);

        auto psidag = dag(psi1);

        psidag.prime("Link");

        auto li_1 = leftLinkIndex(psi1, m);

        auto C = prime(psi1(m), li_1) * op_m;
        C *= prime(psidag(m), "Site");

        for (int k = m + 1; k < j; ++k)
        {
            C *= psi1(k);
            C *= psidag(k);
        }

        auto lj = rightLinkIndex(psi1, j);

        C *= prime(psi1(j), lj) * op_j;
        C *= prime(psidag(j), "Site");

    // Debugging: Print the indices of C to ensure correct contraction

// std::cout << "Indices of C before elt: " << inds(C) << std::endl;

        auto result = elt(C); 
        println("Correlation (<S^z_M S^z_", j, ">) = ", result);
    }
    else if (m > j)
    {
        psi1.position(j);

        auto psidag = dag(psi1);

        psidag.prime("Link");

        auto lj_1 = leftLinkIndex(psi1, j);

        auto C = prime(psi1(j), lj_1) * op_j;
        C *= prime(psidag(j), "Site");

        for (int k = j + 1; k < m; ++k)
        {
            C *= psi1(k);
            C *= psidag(k);
        }

        auto li = rightLinkIndex(psi1, m);

        C *= prime(psi1(m), li) * op_m;
        C *= prime(psidag(m), "Site");

     // Debugging: Print the indices of C to ensure correct contraction

// std::cout << "Indices of C before elt: " << inds(C) << std::endl;

        auto result = elt(C); 
        println("Correlation (<S_z(M) S_z", (j), ">) = ", result);
          }
       }
    }
 }

return 0;
}

After the TDVP time evolution, the state \psi (represented as an MPS) is being printed as

psi =

ITensor ord=2: (dim=2|id=261|“Link”) (dim=2|id=645|“n=1,Site,S=1/2”)
{norm=1.00 (Dense Real)}

ITensor ord=3: (dim=4|id=239) (dim=2|id=65|“n=2,Site,S=1/2”) (dim=2|id=261|“Link”)
{norm=1.41 (Dense Real)}

ITensor ord=3: (dim=8|id=210|“Link”) (dim=2|id=426|“n=3,Site,S=1/2”) (dim=4|id=239)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=193) (dim=2|id=232|“n=4,Site,S=1/2”) (dim=8|id=210|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=377|“Link”) (dim=2|id=282|“n=5,Site,S=1/2”) (dim=4|id=193)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=900) (dim=2|id=752|“n=6,Site,S=1/2”) (dim=8|id=377|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=807|“Link”) (dim=2|id=991|“n=7,Site,S=1/2”) (dim=4|id=900)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=189) (dim=2|id=788|“n=8,Site,S=1/2”) (dim=8|id=807|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=286|“Link”) (dim=2|id=577|“n=9,Site,S=1/2”) (dim=4|id=189)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=667) (dim=2|id=719|“n=10,Site,S=1/2”) (dim=8|id=286|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=356|“Link”) (dim=2|id=469|“n=11,Site,S=1/2”) (dim=4|id=667)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=808) (dim=2|id=160|“n=12,Site,S=1/2”) (dim=8|id=356|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=886|“Link”) (dim=2|id=157|“n=13,Site,S=1/2”) (dim=4|id=808)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=331) (dim=2|id=98|“n=14,Site,S=1/2”) (dim=8|id=886|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=12|“Link”) (dim=2|id=447|“n=15,Site,S=1/2”) (dim=4|id=331)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=951) (dim=2|id=785|“n=16,Site,S=1/2”) (dim=8|id=12|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=8|id=719|“Link”) (dim=2|id=110|“n=17,Site,S=1/2”) (dim=4|id=951)
{norm=2.00 (Dense Real)}

ITensor ord=3: (dim=4|id=360) (dim=2|id=360|“n=18,Site,S=1/2”) (dim=8|id=719|“Link”)
{norm=2.83 (Dense Real)}

ITensor ord=3: (dim=2|id=890|“Link”) (dim=2|id=401|“n=19,Site,S=1/2”) (dim=4|id=360)
{norm=2.00 (Dense Real)}

ITensor ord=2: (dim=2|id=68|“n=20,Site,S=1/2”) (dim=2|id=890|“Link”)
{norm=1.41 (Dense Real)}

I verified the energy against the Exact Diagonalization (ED) results, and they match well. However, when calculating the S_z(M)S_z(i)> correlation function, where M is the middle spin (used as the reference) and i is any other site index, it shows the following error:

Correlation (<S^z_M S^z_>) = 9 1
Correlation (<S^z_M S^z_1>) = 3.78032e-26
Correlation (<S^z_M S^z_>) = 9 3
inds() = (dim=4|id=405)’ (dim=4|id=405)
From line 89, file itensor.cc

Wrong number of IndexVals passed to elt/eltC (expected 0, got 2)

Wrong number of IndexVals passed to elt/eltC (expected 0, got 2)
Aborted (core dumped)

Any suggestion would be really helpful.

Thank you in advance.

miles · December 4, 2024, 8:06pm

It seems your code has a bug in it. You can print out various objects and intermediate steps to find this bug, since many objects in ITensor are printable.

This particular error about “Wrong number of IndexVals” typically means that you thought you code is resulting in a scalar, but that the final tensor you are getting is not a scalar and still has uncontracted indices. So most likely the bug is that you are not contracting indices properly (either contracting some too early, so they disappear and/or not taking account of prime levels correctly).

For your specific calculation, you may also be able to use our correlation_matrix function, as documented on this page:
http://itensor.org/docs.cgi?vers=cppv3&page=classes/mps_mpo_algs

Finally, we’d appreciate it if you format your code with three backticks (```) when posting code

miles · December 5, 2024, 2:22am

Hi Jaga,
Really glad that you were able to figure out the problem and that my comments helped. Also yes, please check out the correlation_matrix function which is quite convenient.

Miles

Sambu_1234 · December 5, 2024, 2:30am

Hi Miles,

Thank you so much for the helpful clues! I was able to find the solution, and the code is now working fine. The issue was related to uncontracted indices, caused by missing links between a few site tensors while creating the initial entangled pairs.

I also truly appreciate the information about the correlation_matrix function and will ensure to format the code with backticks (```) when posting moving forward.

Best,
Sambu

Sambu_1234 · December 5, 2024, 11:46am

Thanks @miles! I used the correlation_matrix function for my calculation, and it is very helpful and efficient, saving time on tasks that would otherwise take much longer.

Sambu

Topic		Replies	Views
Dynamical correlation function for thermal states ITensor Julia Questions	4	308	October 28, 2023
How to derive thermal state by imaginary time evolution (purification) ITensor Julia Questions julia , mps , mpo	3	344	June 19, 2023
Calculating finite temperature two time correlation function for fixed magnetisation sector ITensor Julia Questions	8	396	September 5, 2023
Measurement of local obervables ITensor C++ Questions	5	199	August 21, 2023
Complex values of time dependent single site spin correlations ITensor Julia Questions julia , mps	5	215	October 5, 2023

Sz(M)Sz(i) Correlation Function in a Thermal State

Related topics