Hi miles, thank you again for the help! I am using Julia version 1.7.1, and ITensors v0.2.12

Yes, the DMRG calculation is finished. And then I use my own function to do some correlator measurement. Here I do not use the built-in correlation_matrix because my site type is customized, and correlation_matrix can not work

In case you want to check my measurement function. I attach my code below:

sites = [

if mod(n,3) == 2

Index([QN()=>2]; tags="S=1/2, Site, n=$n")

else

siteind("Fermion", n; conserve_qns=true)

end

for n = 1:3*ns

]

state = [

if mod(n,3) == 2

"Up"

else

(mod(n,fl*3) == 1 || mod(n,fl*3) == 0) ? "Occ" : "Emp"

end

for n = 1:3*ns

]

psii = randomMPS(sites, state)

# Sz-Sz

function szz(psi::MPS)

szz = zeros(ns,ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) == 2

i1 += 1

j1 = i1

orthogonalize!(psi, i)

szz[i1,i1] = scalar(ts*psi[i]*op(szo,sites,i)*prime(op(szo,sites,i))*dag(prime(prime(psi[i],"Site"),"Site")))

lind = linkind(psi,i)

tsi = ts*psi[i]*op(szo,sites,i)*dag(prime(prime(psi[i],"Site"),lind))

for j = i+1:3*ns

tsi *= psi[j]

if mod(j,3) == 2 # spin site only

j1 += 1

lind = linkind(psi,j-1)

szz[i1,j1] = scalar(tsi*op(szo,sites,j)*dag(prime(prime(psi[j],"Site"),lind)))

szz[j1,i1] = conj(szz[i1,j1])

end

tsi *= dag(prime(psi[j],"Link")) # no JW string

end

end

end

return szz

end

# Sz

function sz(psi::MPS)

sz = zeros(ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) == 2

i1 += 1

orthogonalize!(psi,i)

sz[i1] = scalar(ts*psi[i]*op(szo,sites,i)*dag(prime(psi[i],"Site")))

end

end

return sz

end

# Nf

function n(psi::MPS)

n = zeros(2*ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) != 2

i1 += 1

orthogonalize!(psi, i)

n[i1] = scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))

end

end

return n

end

# Nf-Nf

function nn(psi::MPS)

nn = zeros(2*ns,2*ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) != 2

i1 += 1

j1 = i1

orthogonalize!(psi,i)

nn[i1,i1] = scalar(ts*psi[i]*op("N",sites,i)*prime(op("N",sites,i))*dag(prime(prime(psi[i],"Site"),"Site")))

lind = linkind(psi,i)

tsi = ts*psi[i]*op("N",sites,i)*dag(prime(prime(psi[i],"Site"),lind))

for j = i+1:3*ns

tsi *= psi[j]

if mod(j,3) != 2

j1 += 1

lind = linkind(psi,j-1)

nn[i1,j1] = scalar(tsi*op("N",sites,j)*dag(prime(prime(psi[j],"Site"),lind)))

nn[j1,i1] = conj(nn[i1,j1])

tsi *= dag(prime(psi[j],"Link"))

else

tsi *= dag(prime(psi[j],"Link"))

end

end

end

end

return nn

end

# Cdag-C

function ff(psi::MPS)

ff = zeros(2*ns,2*ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) != 2

i1 += 1

j1 = i1

orthogonalize!(psi,i)

ff[i1,i1] = scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))

lind = linkind(psi,i)

tsi = ts*psi[i]*op("Cdag",sites,i)*dag(prime(prime(psi[i],"Site"),lind))

for j in i+1:3*ns

tsi *= psi[j]

if mod(j,3) != 2

j1 += 1

lind = linkind(psi, j-1)

ff[i1,j1] = scalar(tsi*op("C",sites,j)*dag(prime(prime(psi[j],"Site"),lind)))

ff[j1,i1] = conj(ff[i1,j1])

# with JW string

tsi *= op("F",sites,j)*dag(prime(psi[j]))

else

# no JW string

tsi *= dag(prime(psi[j],"Link"))

end

end

end

end

return ff

end

# Fz

# fzi = Nli - Nri / 2

function fz(psi::MPS)

fz = zeros(ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) == 1

i1 += 1

orthogonalize!(psi,i)

nup = scalar(ts*psi[i]*op("N",sites,i)*dag(prime(psi[i],"Site")))

orthogonalize!(psi,i+2)

ndn = scalar(ts*psi[i+2]*op("N",sites,i+2)*dag(prime(psi[i+2],"Site")))

fz[i1] = (nup-ndn)/2

end

end

return fz

end

# Fx & Fy

# sxi = fli fri + fri fli / 2

# syi = -i fli fri + i fri fli / 2

function fxfy(psi::MPS)

fx = zeros(ns)

fy = zeros(ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) == 1

i1 += 1

orthogonalize!(psi,i)

lind1 = linkind(psi,i)

lind2 = linkind(psi,i+1)

tsi = ts*psi[i]*op("Cdag",sites,i)*dag(prime(prime(psi[i],"Site"),lind1))

tsi *= psi[i+1]*dag(prime(psi[i+1],"Link"))

tsi *= psi[i+2]*op("C",sites,i+2)*dag(prime(prime(psi[i+2],"Site"),lind2))

lri = scalar(tsi)

tsi = ts*psi[i]*op("C",sites,i)*dag(prime(prime(psi[i],"Site"),lind1))

tsi *= psi[i+1]*dag(prime(psi[i+1],"Link"))

tsi *= psi[i+2]*op("Cdag",sites,i+2)*dag(prime(prime(psi[i+2],"Site"),lind2))

rli = scalar(tsi)

fy[i1] = (-lri+rli)/2

fx[i1] = (lri+rli)/2

end

end

return fx,fy

end

# Cdga-Sz-C

function rungTunneling(psi::MPS)

rungTunneling = zeros(ns)

ts = ITensor(1.)

i1 = 0

for i in 1:3*ns

if mod(i,3) == 1

i1 += 1

orthogonalize!(psi,i)

lind1 = linkind(psi,i)

lind2 = linkind(psi,i+1)

tsi = ts*psi[i]*op("Cdag",sites,i)*dag(prime(prime(psi[i],"Site"),lind1))

tsi *= psi[i+1]*op(szo,sites,i+1)*dag(prime(psi[i+1]))

tsi *= psi[i+2]*op("C",sites,i+2)*dag(prime(prime(psi[i+2],"Site"),lind2))

rungTunneling[i1] = scalar(tsi)

end

end

return rungTunneling

end