Define `MPOHamiltonian` with open boundary condition

[4psireal2]

I'm not sure how to convert the MPO representation as Vector of the Hamiltonian for ex. the transverse Ising chain into MPSKit's MPOHamiltonian. I think it can be casted into a periodic array but I'm not sure about the details. Are there perhaps examples? Many thanks!

[lkdvos]

There is currently a bit of a weird quirk, in the sense that we are using the same struct for both infinite and finite hamiltonians. It is at the time of the call to find_groundstate etc that the Hamiltonian will know whether or not it needs to be infinite or finite, and then just automatically adapt. With that in mind, you should be able to just input the MPO representation of your hamiltonian into MPOHamiltonian(), however it expects it to be flattened, so a 3-dimensional array where O[i,j,k] is the entry at site i, on row j and column k. I do apologize for the lack of documentation, and will try to resolve this as soon as possible. In the meantime, I would also suggest having a look at MPSKitModels.jl, which has a lot of functionality for easily constructing mpo representations.

[4psireal2]

Thank you for the prompt response! Just to double check, here is my Hamiltonian for the transverse Ising model for 3 sites:

T = ComplexF64
X = TensorMap(T[0 1;1 0], ℂ^2 ← ℂ^2)
X_T = TensorMap(T[0 1;1 0], ℂ^2 ← ℂ^2)
Z = TensorMap(T[1 0;0 -1], ℂ^2 ← ℂ^2)
Id_mat = TensorMap(T[1 0;0 1], ℂ^2 ← ℂ^2)
J = 1.0
g = 1.0

## 3-site Hamiltonian
data = Array{Any,3}(missing,3,3,3)
data[1, 1, 1] = Id_mat
data[1, 1, 2] = -J * X
data[1, 1, 3] = g * Z
data[2, 1, 1]= Id_mat
data[2, 1, 2] = -J * X
data[2, 1, 3] = g * Z
data[2, 2, 3] = X
data[2, 3, 3] = Id_mat
data[3, 1, 1] = g * Z
data[3, 2, 1] = -J * X
data[3, 3, 1] = Id_mat

Is this correct? I then ran the DMRG method:

L = 3
D = 128
init_state = FiniteMPS(L, ℂ^2, ℂ^D)
OBC_Ising = MPOHamiltonian(data)
ψ, envs, δ = find_groundstate(init_state, OBC_Ising, DMRG())

I noticed that the structure of init_state

┌ CL[4]: TensorMap(ℂ^1 ← ℂ^1)
├── AL[3]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^1)
├── AL[2]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^2)
└── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^2)) 

differs from that of ψ

┌── AC[3]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^1)
├── AL[2]: TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^2)
└── AL[1]: TensorMap((ℂ^1 ⊗ ℂ^2) ← ℂ^2)

Why is this the case?

[lkdvos]

Your hamiltonian looks correct, and the different structures you get is actually just a matter of printing. The FiniteMPS actually only lazily keeps all of the AL, AR and AC tensors, in the sense that it only computes them when queried, and then stores that result. When you print a state, it tries to compute the least amount of additional data, which can change depending on what information is available.

[4psireal2]

I obtained very strange energy values for the ground state of the transverse Ising model. Can you perhaps take a look?

""" Build a tranverse field Ising Hamiltonian
Hamiltonian: -J(∑_{i,j} σ^x_i  σ^x_j + g∑_i σ^z_i)   
Guide: https://github.com/maartenvd/MPSKit.jl/blob/45cab94fb3d78169a5987c54ed2b762652dd4c70/docs/src/man/operators.md """

using MPSKit, MPSKitModels, TensorKit
using ProgressMeter, Plots
using LinearAlgebra

L = 16
D = 128

T = ComplexF64
X = TensorMap(T[0 1;1 0], ℂ^2 ← ℂ^2)
Z = TensorMap(T[1 0;0 -1], ℂ^2 ← ℂ^2)
Id_mat = TensorMap(T[1 0;0 1], ℂ^2 ← ℂ^2)
J = 1.0
g = 1.0

init_state = FiniteMPS(L, ℂ^2, ℂ^D)

## L-site Hamiltonian
data = Array{Any,3}(missing,L,3,3)
data[1, 1, 1] = Id_mat
data[1, 1, 2] = X
data[1, 1, 3] = -J*g*Z
for i = 2:L-1
    data[i, 1, 1]= Id_mat
    data[i, 1, 2] = -J * X
    data[i, 1, 3] = -J * g * Z
    data[i, 2, 3] = -J * X
    data[i, 3, 3] = Id_mat
end
data[L, 1, 1] = -J*g*Z
data[L, 2, 1] = X
data[L, 3, 1] = Id_mat

## OBC
println("Computation 1")
OBC_Ising = @mpoham sum(-J * S_xx(){i, j} -J * g * S_z(){k} for (i,j) in nearest_neighbours(FiniteChain(L)) for k in vertices(FiniteChain(L)))
ψ_obc, env_obc, δ_obc = find_groundstate(init_state, OBC_Ising; verbose=true) #XXX: strange energy

println("Computation 2")
OBC_Ising = MPOHamiltonian(data)
ψ_obc_own, env_obc_own, δ_obc_own = find_groundstate(init_state, OBC_Ising; verbose=true) #XXX: strange energy, ϵ = NaN

## 1-site Hamiltonian
data = Array{Any,3}(missing, 1, 3, 3)
data[1, 1, 1] = Id_mat
data[1, 1, 2] = -J * X
data[1, 1, 3] = -J * g * Z
data[1, 2, 3] = -J * X
data[1, 3, 3] = Id_mat

## PBC
println("Computation 3")
PBC_Ising = periodic_boundary_conditions(MPOHamiltonian(data), L)
ψ_pbc_own, env_pbc_own, δ_pbc_own = find_groundstate(init_state, PBC_Ising; verbose=true) #XXX: strange energy, ϵ = NaN

println("Computation 4")
PBC_Ising = periodic_boundary_conditions(transverse_field_ising(; g=g), L)
ψ_pbc, env_pbc, δ_pbc = find_groundstate(init_state, PBC_Ising; verbose=true) #XXX: strange energy

println("Literature energy E/$L: -1.2510242438)") #Ref: https://dmrg101-tutorial.readthedocs.io/en/latest/tfim.html
println("OBC Ising own computation: $δ_obc_own") # δ=2.0e-12
println("OBC Ising MPSKit computation: $δ_obc") # δ=2.0e-12
println("PBC Ising own computation: $δ_pbc_own") # δ=2.0e-12
println("PBC Ising MPSKit computation: $δ_pbc") # δ=2.0e-12

[Gertian]

Hi @4psireal2 ,

I don't have much experience with periodic boundary conditions but when working with open boundary conditions I get this (I included quite a bit of information in my comments in the code below) :

using MPSKit, MPSKitModels, TensorKit
using ProgressMeter, Plots
using LinearAlgebra
using QuadGK

#first a note on the GS energy of the ising model : https://theory.leeds.ac.uk/interaction-distance/applications/ising/map-to-free/
#H = sum_i -σ^x_i σ^x_{i+1} + hσ^z_i which in the thermodynamic limit can be diagonalized by a Jordan-Wigner transformation to a quadratic fermionic Hamiltonian.
#The dispersion relation of this Hamiltonian is E(k) = ±sqrt(1 + h^2 - 2hcos(k)) which is critical when g = 1.
#At h=1 the GS energy is -int_0^pi dk E(k) which is :
h = 1 #critical point
E_GS = -quadgk(k -> sqrt(1 + h^2 - 2*h*cos(k)), -pi, pi)[1]/(2*pi)

#Let us define a finite lattice
L = 100
lattice = FiniteChain(L)    
#Let us define an initial state with bond dimension D on this finite lattice :
D = 10
init_state = FiniteMPS(L, ℂ^2, ℂ^D)

#1) use @mpoham to define the Hamiltonian on this finite lattice
NN_part = @mpoham sum(-4*S_xx(){i,j} for (i,j) in nearest_neighbours(lattice)    )    #note the 4*S_xx ! This is because S_xx is the spin matrices but if you want the critical point of Ising to appear at h = J = 1 we need to work with Pauli matrices !
on_site = @mpoham sum(+h*2*S_z(){i} for i in vertices(lattice)) #same comment !
OBC_Ising = NN_part + on_site #add the two MPO hamiltonians. Surely you could have used a one-liner to make this but I'm not sure it would have been clearer...

#2) manually construct the Hamiltonian
#First we will need the Pauli matrices 
T = ComplexF64
X = TensorMap(T[0 1;1 0], ℂ^2 ← ℂ^2)
Z = TensorMap(T[1 0;0 -1], ℂ^2 ← ℂ^2)

##The finite state machine representation of the Ising Hamiltonian is :
data = Array{Any,3}(missing,L,3,3)
data[1, :, :] = [1 -X h*Z ; missing missing missing ; missing missing missing] 
for i = 2:L-1
    data[i, :, :]= [1 -X h*Z ; missing missing X ; missing missing 1]
end
data[L, :, :] = [missing missing h*Z ; missing missing missing ; missing missing 1]
#from this we get : 
OBC_Ising_manual = MPOHamiltonian(data)

#Let us now calculate the GSs for these things :
ψ, env, δ = find_groundstate(init_state, OBC_Ising; verbose=true) #ψ, env, δ are respectively the groundstate, the environment and the convergence error of the groundstate search
ψ_manual, env_manual, δ_manual = find_groundstate(init_state, OBC_Ising_manual; verbose=true) #ψ, env, δ are respectively the groundstate, the environment and the convergence error of the groundstate search

#and finally, let us compare the aquired energies with those of the analytical solution :
E_density_ising = expectation_value(ψ, OBC_Ising, env) 
@show E_density_ising.-E_GS #All right this matches well in the bulk ! Near the edges we see some deviation which is to be expected on a finite lattice with OBC...

E_density_ising_manual = expectation_value(ψ_manual, OBC_Ising_manual, env_manual)
@show E_density_ising_manual.-E_GS #All right this matches well in the bulk ! Near the edges we see some deviation which is to be expected on a finite lattice with OBC...

Your main mistake was that you used the spin matrices instead of the Pauli matrices in your definition of the Hamiltonian (and were therefor not calculating the energy density at criticality) . Apart from this it seems that you were interpreting the third output of find_groundstate as the final energy which is not true (it's the error on the convergence of the algorithm). Instead you should use : expectation_value(state, Hamiltonian, environments) to get a list of all energy expectation values.

I hope this helps :)

[4psireal2]

Many thanks for the helpful answer @Gertian ! Another related question: I want to define a definite state

init_state_MPS = FiniteMPS(L, ℂ^2, ℂ^χ) 

basis_0 = TensorMap([1, 0], ℂ^2 ← ℂ^1)
basis_1 = TensorMap([0, 1], ℂ^2 ← ℂ^1)

init_product_state = basis_0 ⊗ basis_0 ⊗basis_0 ⊗ basis_1 ⊗ basis_1 ⊗ basis_1 
# TensorMap((ℂ^2 ⊗ ℂ^2 ⊗ ℂ^2 ⊗ ℂ^2 ⊗ ℂ^2 ⊗ ℂ^2) ← (ℂ^1 ⊗ ℂ^1 ⊗ ℂ^1 ⊗ ℂ^1 ⊗ ℂ^1 ⊗ ℂ^1))

Do I now need to perform a SVD so that init_product_state has a similar MPS representation as the init_state_MPS?

Also I'm aware that this won't make a difference but more accurately, it should be as follows right?

for i = 2:L-1
    data[i, :, :]= [1 -X h*Z ; missing missing transpose(X) ; missing missing 1]
end

P/S: I apologize if this is not the place for asking specific questions.

[Gertian]

@4psireal2 I don't see how you could do this with a singular value decomposition since your init_product_state is a map from ℂ^6 to ℂ^2^6 whereas the desired quantum state is a vector on ℂ^2^6 .

I'm guessing you wanted the ℂ^1 legs to represent the virtual dimension but then you should remember that an MPS tensor always has 3 legs such that we need two such legs ! With this in mind we get :

basis_0 = TensorMap([1, 0. *im], ℂ^1 ⊗  ℂ^2 ← ℂ^1) 
basis_1 = TensorMap([0. *im, 1], ℂ^1 ⊗  ℂ^2 ← ℂ^1)
FiniteMPS( [basis_0, basis_0, basis_0, basis_1, basis_1, basis_1]; normalize = true )

Note my somewhat weird notation [1, 0. *im] , this ensures that our state ends up being a vector with complex float values rather than real integer values.

As for your comment on transpose(X), this is indeed entirely irrelevant when you're working with tensormaps over ℂ but when you would implement symmetries this would become relevant. I'm not sure if simply transposing would be enough in this case though.

Honestly I also don't know if this is the place to ask questions like this. @lkdvos I'm guessing it's useful since other people might later get the same problems and use these discussions as a guideline ?

[lkdvos]

Honestly, I see no problems with having these discussions here, it's also helpful for me to remember to make the documentation clearer.

That being said, as far as I understand, you could indeed define the basis states like you did, but actually you would want them to be Tensor instead of TensorMap, which is equivalent with replacing rhe C^1 with C^0. This is the generalized version of a vector, and then you can take their tensor product like you did, which would create the sense state vector. There is a function that then turns this into an MPS, but I think it's not exported: MPSKit.decompose_local_mps will take a tensor and spit out a vector of mpstensors which you can then feed to the finitemps constructor. However, because of the whole mps-tensor needs 3 legs thing, it will interpret the first and last space as utility legs of the MPS, so you would probably want a trivial space there. I'm currently out of office, but I'll try and add a code snippet when I get back. This method is however not very scaleable, as it does construct the dense vector first, and then decomposes it with SVDs. in that sense, Gertian's solution is a bit better.

[Jutho]

Small correction to @lkdvos' comment about the Tensor vs TensorMap issue: this is not achieved by ℂ^0, as this is really a zero-dimensional space, and leads to total dimension zero. ℂ^2 <- ℂ^0 is like a matrix with two rows and zero columns. It does not have any entries.

Instead, what you intended is Tensor([1, 0], ℂ^2), which is equivalent to TensorMap([1, 0], ℂ^2 ← (ℂ^1)^0). The idea is that a tensor product of zero factors of vector spaces is like having zero additional indices. Indeed, you can check that dim(ℂ^0) == 0 whereas dim((ℂ^1)^0) == 1, and the same for any (ℂ^k)^0 for that matter. However,(ℂ^k)^0 is not really the recommend syntax for constructing a tensor product of zero factors. Since the tensor product of vector spaces is also achieved with the ordinary multiplication symbol between vector spaces for convenience, the tensor product of zero factors corresponds to the "multiplicative identity" and can be obtained as one(ℂ^k).