Large spin DMRG

[Lee]

Hi everyone!

I’d like to ask a question about Julia 1.11. This version only supports spin-1 and spin-1/2 calculations and doesn’t work for antiferromagnetic one-dimensional chains with large spins.
I’ve tried the custom definition method shown below, but I still can’t run the ground state calculation correctly using DMRG. Does anyone know what might be causing this issue?
Any suggestions would be greatly appreciated!

using ITensors, ITensorMPS, Statistics, Serialization

import ITensors: siteinds, siteind, op  function siteind(::SiteType"Spin3Half", n::Int; kwargs...)
    return Index(4; tags="Site,Spin3Half,n=$n", kwargs...)
end

function siteinds(::Val{:Spin3Half}, N::Int; kwargs...)
    return [siteind(SiteType("Spin3Half"), i; kwargs...) for i in 1:N]
end

function op(::OpName"Sz", ::SiteType"Spin3Half", s::Index)
    sz = ITensor(prime(s), s)
    sz[prime(s)[1], s[1]] = 3/2
    sz[prime(s)[2], s[2]] = 1/2
    sz[prime(s)[3], s[3]] = -1/2
    sz[prime(s)[4], s[4]] = -3/2
    return sz
end

function op(::OpName"S+", ::SiteType"Spin3Half", s::Index)
    S = 3/2
    sp = ITensor(prime(s), s)
    
    sp[prime(s)[2], s[1]] = sqrt(S*(S+1) - (1/2)*(3/2))  # √3（|1/2⟩→|3/2⟩）
    sp[prime(s)[3], s[2]] = sqrt(S*(S+1) - (-1/2)*(1/2)) # 2.0（|-1/2⟩→|1/2⟩）
    sp[prime(s)[4], s[3]] = sqrt(S*(S+1) - (-3/2)*(-1/2))# √3（|-3/2⟩→|-1/2⟩）
    return sp
end

function op(::OpName"S-", ::SiteType"Spin3Half", s::Index)
    return dag(op(OpName("S+"), SiteType("Spin3Half"), s))
end

function op(::OpName"Sx", ::SiteType"Spin3Half", s::Index)
    sp = op(OpName("S+"), SiteType("Spin3Half"), s)
    sm = op(OpName("S-"), SiteType("Spin3Half"), s)
    return (sp + sm) / 2
end

function op(::OpName"Sy", ::SiteType"Spin3Half", s::Index)
    sp = op(OpName("S+"), SiteType("Spin3Half"), s)
    sm = op(OpName("S-"), SiteType("Spin3Half"), s)
    return (sp - sm) / (2im)
end

let
    N = 72                      
    J = 1.0              
    mps_filename = "spin3half_groundstate.mps"  

    sites = siteinds("Spin3Half", N)

    os = OpSum()
    for j in 1:N
        j_next = j < N ? j+1 : 1
        γ = j % 3 == 1 ? "Sx" : (j % 3 == 2 ? "Sz" : "Sy")
        os += J, "Sx", j, "Sx", j_next
        os += J, "Sy", j, "Sy", j_next
        os += J, "Sz", j, "Sz", j_next
    end
    H = MPO(os, sites)

[miles]

When you say you can’t run the ground state calculation correctly, do you mean that it causes an error to happen? Or that you get a result that you believe is incorrect? Thanks.

[Lee]

Dear Dr. Miles,

Thank you very much for your help.

To answer directly: the problem is that the calculated result appears incorrect to me, rather than an explicit error during runtime. When simulating spin-1/2 chain, my results align well with the theoretical predictions from this work and similar literature. However, for spin-3/2 chain, the code yields a ground state energy of approximately -2.1, which deviates significantly from the theoretical value of -2.8.For context, I referenced the study Conformal anomaly and critical exponents of Heisenberg spin models with half-integer spin (DOI: https://doi.org/10.1103/PhysRevB.36.8582) for validation.

Given this substantial discrepancy, I suspect there may be an issue with my code implementation, especially considering Julia 1.11’s native support is limited to spin-1 calculations. I would greatly appreciate your insights or suggestions on potential causes of this mismatch.

Thank you again for your time and assistance.

[miles]

Thanks. First off, I want to clear up some possible confusion, which is that the support for certain kind of spin such as spin-1 isn’t related to Julia but to a particular library, ITensorMPS, that is written in the Julia programming language. ITensorMPS defines a S=1 site type by default as you know, and this support is not tied to any particular version of the Julia programming language.

As for the energy mismatch, could it just be coming from finite-size effects or convergence issues? What is the system size used for the -2.8 energy in that paper? Is that an estimate of the infinite-system-size energy (bulk energy)? What system size(s) did you use? With an open boundary system there can be substantial boundary contributions to the energy and it can be imporatant to carefully extrapolate these away or subtract them off.

[Lee]

Dear Dr. Miles,

Thank you very much for your prompt reply and for clarifying the relationship between ITensorMPS and Julia in terms of spin support—your explanation has helped me understand this clearly.

Regarding the energy mismatch issue, I would like to share two additional observations after further verification: First, for spin-1/2 calculations, even with small system sizes, the energy interpolation does not show such a significant deviation. Second, for spin-3/2, we used a system with 72 sites under periodic boundary conditions and set maxdim to 1400, but the calculated energy still stabilized around 2.1, with a persistent gap from the expected result. Therefore, I still suspect there may be oversights in the dimension definition or calling method of the spin operators in my code.

Given this, would you be kind enough to provide a code example for the extended definition of spin-3/2 or higher-spin operators? It would be incredibly helpful for me to troubleshoot the problem.

Thank you again for your guidance and support amidst your busy schedule. I look forward to your response.

[mtfishman]

I see you’ve defined:

function op(::OpName"S-", ::SiteType"Spin3Half", s::Index)
    return dag(op(OpName("S+"), SiteType("Spin3Half"), s))
end

That doesn’t look correct since S- should be the transpose of S+, while dag only performs complex conjugation and flipping arrow directions (related to a transpose but in ITensor not the same thing, since you also have to swap the row and column index labels). I think this should be:

function op(::OpName"S-", ::SiteType"Spin3Half", s::Index)
    return swapprime(dag(op(OpName("S+"), SiteType("Spin3Half"), s)), 0 => 1)
end

[Lee]

Dear mtfishman,

Thank you so much for pointing out the issue with my S⁻ operator definition in ITensor. Your note about S⁻ needing to be the transpose of S⁺—instead of just using dag, which misses swapping row/column indices—helped me fix the misunderstanding quickly.

To avoid further errors, I explicitly wrote out the matrixs (see code below), and now the Hamiltonian and MPO run correctly without issues.

module SpinFiveHalvesSite
using ITensors, ITensorMPS

ITensors.space(::SiteType"S=5/2") = 6  

ITensors.op(::OpName"Sz", ::SiteType"S=5/2") = [
    2.5  0.0   0.0    0.0    0.0    0.0
    0.0  1.5   0.0    0.0    0.0    0.0
    0.0  0.0   0.5    0.0    0.0    0.0
    0.0  0.0   0.0   -0.5    0.0    0.0
    0.0  0.0   0.0    0.0   -1.5    0.0
    0.0  0.0   0.0    0.0    0.0   -2.5
]

ITensors.op(::OpName"S+", ::SiteType"S=5/2") = [
    0.0  sqrt(5)   0.0    0.0     0.0      0.0
    0.0   0.0    sqrt(8)  0.0     0.0      0.0
    0.0   0.0     0.0     3.0     0.0      0.0
    0.0   0.0     0.0     0.0    sqrt(8)   0.0
    0.0   0.0     0.0     0.0     0.0     sqrt(5)
    0.0   0.0     0.0     0.0     0.0      0.0
]

ITensors.op(::OpName"S-", ::SiteType"S=5/2") = [
    0.0    0.0     0.0     0.0      0.0      0.0
    sqrt(5) 0.0     0.0     0.0      0.0      0.0
    0.0   sqrt(8)   0.0     0.0      0.0      0.0
    0.0    0.0     3.0     0.0      0.0      0.0
    0.0    0.0     0.0    sqrt(8)   0.0      0.0
    0.0    0.0     0.0     0.0     sqrt(5)   0.0
]


ITensors.op(::OpName"Sx", ::SiteType"S=5/2") = [
    0.0       sqrt(5)/2    0.0      0.0       0.0        0.0
    sqrt(5)/2   0.0      sqrt(8)/2  0.0       0.0        0.0
    0.0       sqrt(8)/2   0.0      3.0/2      0.0        0.0
    0.0        0.0       3.0/2     0.0      sqrt(8)/2   0.0
    0.0        0.0        0.0     sqrt(8)/2   0.0      sqrt(5)/2
    0.0        0.0        0.0       0.0      sqrt(5)/2   0.0
]


ITensors.op(::OpName"Sy", ::SiteType"S=5/2") = [
    0.0        -im*sqrt(5)/2   0.0          0.0           0.0            0.0
    im*sqrt(5)/2   0.0        -im*sqrt(8)/2  0.0           0.0            0.0
    0.0        im*sqrt(8)/2    0.0         -im*3.0/2       0.0            0.0
    0.0          0.0          im*3.0/2      0.0         -im*sqrt(8)/2   0.0
    0.0          0.0           0.0        im*sqrt(8)/2    0.0          -im*sqrt(5)/2
    0.0          0.0           0.0           0.0        im*sqrt(5)/2    0.0
]


end 

if abspath(PROGRAM_FILE) == @__FILE__
    using .SpinFiveHalvesSite  
    using ITensorMPS, ITensors

    s = siteind(SiteType("S=5/2"), 1)
    Sz = op("Sz", SiteType("S=5/2"))
    println("Sz operator for S=5/2 (spin 2.5):")
    println(Sz)
end

Thanks again for your valuable help!