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！

0 votes

Hi JR,

Thanks for the question. In Julia this function is called `contract`

. Here is the documentation on it:

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

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

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)

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)

...