# How to get the MPO of L^+L from the one of L ?[julia]

+1 vote

Hi,
I am trying to get the MPO representation of L^+L from the one of L using the
nmultMPO function, which can be simulated in C++ version. When I use the julia version, iI can not find something like that.
Thanks！

Hi JR,
Thanks for the question. In Julia this function is called contract. Here is the documentation on it:

https://itensor.github.io/ITensors.jl/dev/MPSandMPO.html#ITensors.NDTensors.contract-Tuple{MPO,%20MPO}

Please let us know if you don't get the results you expect. A good way to test that this is working as you want is to make random MPS using the randomMPS function and compute matrix elements of your squared MPO as well as in an exact way using the function inner(A,psi,B,psi) where A and B are MPOs. The function inner(A,psi,B,psi) computes <psi|A^+ B|psi> using a different algorithm so would provide a good check.

Best,
Miles

commented by (550 points)
thanks a lot
commented by (550 points)
1 using ITensors
2 using Printf
3 using Random
4
5 Random.seed!(1234)
6
7 let
8   N = 12
9   weigth=10
10   # Create N spin-one degrees of freedom
11   sites = siteinds("S=1",N)
12   # Alternatively can make spin-half sites instead
13   #sites = siteinds("S=1/2",N)
14
15   # Input operator terms which define a Hamiltonian
16   ampo = AutoMPO()
17   for j=1:N-1
18     ampo += "Sz",j,"Sz",j+1
19     ampo += 0.5,"S+",j,"S-",j+1
20     ampo += 0.5,"S-",j,"S+",j+1
21   end
22   # Convert these terms to an MPO tensor network
23   H = MPO(ampo,sites)
24   cutoff1=1E-8
25   maxdim1=100000
26   H1=contract(H,H;cutoff=cutoff1,maxdim=maxdim1)
27   # Create an initial random matrix product state
28 end

I tested the above codes, and it failed.

ERROR: LoadError: MethodError: no method matching IndexSet(::Nothing, ::Nothing)
Closest candidates are:
IndexSet(::Any) at /Users/mac/.julia/packages/ITensors/Ligbx/src/indexset.jl:115
Stacktrace:
[1] eigen(A::ITensor{0}, Linds::IndexSet{0, Union{}, Tuple{}}, Rinds::IndexSet{0, Union{}, Tuple{}}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:ishermitian, :ortho, :tags, :cutoff, :maxdim, :mindim), Tuple{Bool, String, TagSet, Float64, Int64, Int64}}})
@ ITensors ~/.julia/packages/ITensors/Ligbx/src/decomp.jl:206
[2] factorize_eigen(A::ITensor{2}, Linds::Tuple{}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{5, Symbol}, NamedTuple{(:ortho, :tags, :cutoff, :maxdim, :mindim), Tuple{String, TagSet, Float64, Int64, Int64}}})
@ ITensors ~/.julia/packages/ITensors/Ligbx/src/decomp.jl:337
[3] factorize(A::ITensor{2}, Linds::Tuple{}; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{5, Symbol}, NamedTuple{(:ortho, :tags, :cutoff, :maxdim, :mindim), Tuple{String, TagSet, Float64, Int64, Int64}}})
@ ITensors ~/.julia/packages/ITensors/Ligbx/src/decomp.jl:420
[4] *(A::MPO, B::MPO; kwargs::Base.Iterators.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:cutoff, :maxdim), Tuple{Float64, Int64}}})
@ ITensors ~/.julia/packages/ITensors/Ligbx/src/mps/mpo.jl:518
[5] #contract#588
@ ~/.julia/packages/ITensors/Ligbx/src/mps/abstractmps.jl:1158 [inlined]
[6] top-level scope
@ ~/Documents/dmrg_julia/dmrg.jl:26
in expression starting at /Users/mac/Documents/dmrg_julia/dmrg.jl:7
commented by (14.1k points)
You would have to do:

H2 = contract(prime(H), H; cutoff=cutoff1, maxdim=maxdim1)

to get H^2 or:

Hdag = swapprime(dag(H), 0 => 1)
H2 = contract(prime(Hdag), H; cutoff=cutoff1, maxdim=maxdim1)

to get H^dag H (where I'm using H^dag to denote the conjugate transpose of H).

Note that the result of both operations will be an MPO with site indices with prime levels of 0 and 2, so you would probably want to change them back to 0 and 1 with:

H2 = replaceprime(H2, 2 => 1)
commented by (550 points)