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

Question

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.
Who knows the answer?
Thanks！

Answer 1

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

Comments

thanks a lot

---

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

---

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)

---

Thank you for your answers!  Now that the problem is solved.

---

Now that the problem is solved. Thank you very munch !