Expectation value of an operator after dmrg using ground state.

DMRG and Numerical Methods
1

Hello there to the community hope all are doing good,
I have done a DMRG calculation of 1d spin 1/2 Heisenberg model and found its ground state as an mps, now I wish to calculate some expectations like Sx, Sy, Sz, S.S . The problem I am facing is I am getting the ket vector and its bra as a 2x2 matrix,
my code reads like this:-
N=2
sites = siteinds(“S=1/2”,N)
psi = psi0; #the ground state called from dmrg
orthogonalize!(psi,j)
psidag_j = dag(prime(psi[j], “Site”))
println(psi[j])
println(psidag_j)
Szjop = op(sites, “Sz” , j)
println(Szjop)

ITensor ord=2
Dim 1: (dim=2|id=575|“S=1/2,Site,n=1”)
Dim 2: (dim=2|id=467|“Link,l=1”)
NDTensors.Dense{Float64, Vector{Float64}}
2×2
-0.029907032140422538 0.40802914179299793
-0.9123808978501717 -0.013374831375957739
ITensor ord=2
Dim 1: (dim=2|id=575|“S=1/2,Site,n=1”)’
Dim 2: (dim=2|id=467|“Link,l=1”)
NDTensors.Dense{Float64, Vector{Float64}}
2×2
-0.029907032140422538 0.40802914179299793
-0.9123808978501717 -0.013374831375957739
ITensor ord=2
Dim 1: (dim=2|id=884|“S=1/2,Site,n=1”)’
Dim 2: (dim=2|id=884|“S=1/2,Site,n=1”)
NDTensors.Dense{Float64, Vector{Float64}}
2×2
0.5 0.0
0.0 -0.5

The problem begins when I print the expectation value of the operator;
Sz_expect = scalar(psidag_j * op(sites, “Sz”, j) * psi[j])

The output shows this error;

DimensionMismatch: In scalar(T) or T, ITensor T is not a scalar (it has indices ((dim=2|id=575|“S=1/2,Site,n=1”)‘, (dim=2|id=884|“S=1/2,Site,n=1”)’, (dim=2|id=884|“S=1/2,Site,n=1”), (dim=2|id=575|“S=1/2,Site,n=1”))).
Please help as I am unable to resolve my problem

2

Hi @shuklavaibhav, please next time format better your question using markdowns.

Why your code doesn't work

Anyway, the problem is that you are using a different IndexSet for the MPS and the operator S_z, infact the index id of the site of psi and of S_z do not match, you have to use the same site indices, you can get them from psi as follows:

sites = siteinds(psi0)

By changing this line your code should work.

I would recommend you to use the utility function expect (MPS and MPO · ITensors.jl) you can use it as follows:

mag_profile = expect(psi,"Sz")

by default it will compute the whole magnetization profile, but if you are only interested to a particular site, or range, you can use the sites keyword to select them

1 Like
3

Sorry for the inconvenience I was not aware about markdown will keep in mind next time. Thankyou for your help it worked

4

Can you please also tell how to compute the expectation of “Sy” and “Sy[i]Sy[j]” correlations. As it gives an error :-
InexactError: Float64(0.0 + 0.5im)
for the expectation of “Sy”

5

@shuklavaibhav
Again, please use markdowns, you can find instructions on how to use them here (Please Read: Make It Easier to Help You)

Anyway, are you sure you actually have to compute them? What you are attempting to do is compute the mean value of a the matrix Sy over a real state, this is always going to be zero (this is also true for the non diagonal elements of the correlation matrix C_{i,j}=\langle S_i^yS_j^y\rangle, while the diagonal elements would be the identity since (S_i^y)^2 = 1).

Mathematical demonstration

The mean value of S_y in the computational basis is represented as:

\begin{pmatrix} \alpha^* & \beta^* \end{pmatrix} \begin{pmatrix} 0 & -i\\ i & 0 \end{pmatrix} \begin{pmatrix} \alpha \\ \beta \end{pmatrix}

By performing the multiplications you arrive that this is equal to 2 Re(i\beta^* \alpha), if both \alpha and \beta are real, than i\beta^*\alpha is a purely imaginary number, so its real part is zero.

So if you know that your state is real (and you know that because the S=1/2 Heisenberg Hamiltonian is real and symmetric) you already know that this mean value is zero.

If you still want to do it numerically, you have to convert the state to a complex one

psi = convert_leaf_eltype(ComplexF64, psi)

To compute the correlations the function to use is correlation_matrix(psi, "Sy", "Sy"), please refer to the example (MPS and MPO Examples · ITensors.jl) and to the documentation of the function (MPS and MPO · ITensors.jl)

6

Hi there, thankyou for your valuable guidance, can you help me through calculating the expectation value of the spin spin operator, I saw the code for this in C++ but am not able to write it in julia. I have attached the C++ code link below.
Thank you.

\vec{S}_j \cdot \vec{S}_{j+1} = S^z_j S^z_{j+1} + \frac{1}{2} S^+_j S^-_{j+1} + \frac{1}{2} S^-_j S^+_{j+1}
7

It isn’t much different

TotalSdS = let TotalSdS = 0.0
    for b in 1:(N-1)
        orthogonalize!(psi, b)

        bondket = psi[b] * psi[b+1]
        bondbra = dag(prime(bondket; tags="Site"))

        zzop = op(sites, "Sz", b) * op(sites, "Sz", b + 1)
        pmop = 0.5 * op(sites, "S+", b) * op(sites, "S-", b + 1)
        mpop = 0.5 * op(sites, "S-", b) * op(sites, "S+", b + 1)

        zz = scalar(bondbra * zzop * bondket)
        pm = scalar(bondbra * pmop * bondket)
        mp = scalar(bondbra * mpop * bondket)
        TotalSdS += zz + pm + mp
    end
    TotalSdS
end
@show TotalSdS
8

Thank you for replying and sorry for the inconvenience but it again shows the same error that I started with regarding dimension mismatch
DimensionMismatch: In scalar(T) or T[], ITensor T is not a scalar (it has indices ((dim=2|id=614|"S=1/2,Site,n=1")', (dim=2|id=390|"S=1/2,Site,n=2")', (dim=2|id=151|"S=1/2,Site,n=1")', (dim=2|id=151|"S=1/2,Site,n=1"), (dim=2|id=306|"S=1/2,Site,n=2")', (dim=2|id=306|"S=1/2,Site,n=2"), (dim=2|id=614|"S=1/2,Site,n=1"), (dim=2|id=390|"S=1/2,Site,n=2"))).

9

can you please also post the code and the backtrace?

EDIT: I looked better at your code and I think I found the problem even without the code and the backtrace. I think the sites argument you are passing to the op function, is different from the indices used by psi, you already did the same error before. What I think is not clear is that each index has an ID and a prime level, when you multiply 2 tensors, the contracted indices are the ones that have same ID and same prime level. So if you pass different inds to op respect to the ones used by psi the ID will not match and the indices will not be contracted

10

Can you tell me what is the backtrace I’ll send the code asap

11

What you posted is the error message, but with it you also have a lot of other lines tells you which functions have been called and in which libraries, that is the backtrace.

12 
using ITensors, ITensorMPS
N=4
sites = siteinds("S=1/2",N)
os= OpSum()

for  j in 1:N-1
    os += 0.25, "S+", j, "S-", j+1
    os += 0.25, "S-", j, "S+", j+1
    os += 0.5,"Sz", j, "Sz", j+1
end    

H = MPO(os,sites)

nsweeps = 10 
maxdim = [10,20,50,100,105] # gradually increase states kept 105
cutoff = [1E-10] # desired truncation error
noise = [1E-6]
weight = 20


psi0_init = randomMPS(sites;linkdims=2)
energy0,psi0 = dmrg(H,psi0_init;nsweeps,maxdim,cutoff,noise)

println()

psi=psi0
TotalSds = let TotalSds = 0.0
    for j in 1:N-1
        orthogonalize!(psi, j)
        
        bondket = psi[j] * psi[j+1]
        bondbra = dag(prime(bondket; tags="Site"))
            
        zzop = op(sites, "Sz" , j) * op(sites, "Sz" , j + 1)
        pmop = 0.5 * op(sites, "S+" , j) * op(sites, "S-", j + 1)
        mpop = 0.5 * op(sites, "S-", j) * op(sites, "S+", j + 1)
        
        zz = scalar(bondbra * zzop * bondket)
        pm = scalar(bondbra * pmop * bondket)
        mp = scalar(bondbra * mpop * bondket)
        TotalSds += zz + pm + mp
    end
    TotalSds
end
@show TotalSds
13 
Stacktrace:
 [1] getindex
   @ C:\Users\viainfo\.julia\packages\ITensors\FpnkY\src\itensor.jl:1123 [inlined]
 [2] scalar(T::ITensor)
   @ ITensors C:\Users\viainfo\.julia\packages\ITensors\FpnkY\src\itensor.jl:1030
 [3] top-level scope
   @ .\In[183]:13
14

I somewhat understood what you are saying but what must I do to resolve this?

15

Is this the code you are currently using? cause I copy pasted it and had an error cause os was being used in a local scope while a global variable with the same name was present, so I had to use a let block. After changing this I didn’t get any errors

os = let os = OpSum()
    for j in 1:N-1
        os += 0.25, "S+", j, "S-", j + 1
        os += 0.25, "S-", j, "S+", j + 1
        os += 0.5, "Sz", j, "Sz", j + 1
    end
    os
end
16

Thank you that did work…