Computing Central Charge from Entanglement Entropy

Hi, I’m new to Julia and ITensor so please forgive me if I’m missing a very obvious error in my code.

I’m trying to compute the central charge of the tricritical Ising model at a critical point from the entanglement entropy using the following formula:

S_A = (c/3)\log{((L/\pi)\sin{(\pi l / L)}})

where c is the central charge, l is the subsystem size, and L is the system size. This is (1.1) from Pasquale Calabrese and John Cardy J. Stat. Mech. (2004)

Following “Confusion about the central charge extracted from Entanglement Entropy”, I wrote the following code:

using ITensors
using Plots

function chrd_dist(j, L)
    return L/π * sin(π/L * j)

function entropy(psi)
    psic = copy(psi)
    N = length(psi)
    ent = Vector{Float64}(undef, N-1)
    for jj in 1:N-1
        orthogonalize!(psic, jj)
        U, S, V = svd(psic[jj], (linkind(psic, jj - 1), siteind(psic, jj)))
        ent[jj] = 0.0
        for n=1:dim(S,1)
            p = S[n,n]^2
            ent[jj] -= p * log(p)
    return ent, map(b -> chrd_dist(b, N), 1:N-1)

g = 0.856
N = 100
sites = siteinds("S=1/2", N)

weight = abs(40 * g)

os = OpSum()
for j in 1:N
    os += -2*2, "Sx", j
for j in 1:N-1
    os += -2*4, "Sz", j, "Sz", j+1
for j in 1:N-2
    os += 8*g, "Sx", j, "Sz", j+1, "Sz", j+2
for j in 1:N-2
    os += 8*g, "Sz", j, "Sz", j+1, "Sx", j+2
os += -2*4, "Sz", N, "Sz", 1 # Periodic BC
os += 8*g, "Sx", N-1, "Sz", N, "Sz", 1
os += 8*g, "Sx", N, "Sz", 1, "Sz", 2
os += 8*g, "Sz", N-1, "Sz", N, "Sx", 1
os += 8*g, "Sz", N, "Sz", 1, "Sx", 2

# Specify the site indices explicitly to avoid indexing issues
H = MPO(os, sites, 1:N)

nsweeps = 5
maxdim = [10, 20, 100, 100, 200]
cutoff = [1E-10]

psi0 = randomMPS(sites, 10)

energy, psiof = dmrg(H, psi0; nsweeps, maxdim, cutoff)

entf, chrdf = entropy(psiof)

pl = plot(log.(chrdf), log.(entf), label=false, seriestype=:scatter)
alg = (x, δ, β) -> δ * x + β
c = 1/2
plot!(pl, x -> alg(x, c/3, -0.6), log.(chrdf), label="c = $c", legend=:topleft)

xlabel!(pl, "Log(d(l))")
ylabel!(pl, "S(l)")
title!("Central Charge")


But I keep getting the error:

MethodError: no method matching _indices(::Tuple{}, ::Nothing)

When I do not have a for loop that sums over all the possible bonds, I get no error. Any idea how to fix this?


At site jj = 1, there is no link index linkind(psic, jj - 1) since you are the edge of the system, so the line:

        U, S, V = svd(psic[jj], (linkind(psic, jj - 1), siteind(psic, jj)))

won’t work. You could use U, S, V = svd(psic[jj], (siteind(psic, jj),)) when jj == 1.

Great, thank you very much. I thought this was the problem initially and tried to address it with an if statement but failed… It seems to work now though!

An alternative to using an if-statement would be to use:

U, S, V = svd(psic[jj], (linkinds(psic, jj - 1)..., siteinds(psic, jj)...))

which works because linkinds(psic, 0) outputs an empty collection of indices. Also it is more general since it handles cases such as if certain MPS tensors have 0 or more than 1 site index, or 0 or more than 1 link index on other links besides on the edge of the system.