# MPO construction from two site custom op

**URL:** <https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924>\
**Category:** ITensor Julia Questions\
**Tags:** mpo\
**Created:** [August 13, 2024, 6:41pm UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924 "2024-08-13T18:41:38Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bhargava\_Balaganchi](https://avatars.discourse-cdn.com/v4/letter/b/a5b964/32.png) [@Bhargava\_Balaganchi](https://itensor.discourse.group/u/Bhargava_Balaganchi)\
**Post date:** [August 13, 2024, 6:41pm UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/1 "2024-08-13T18:41:38Z")

</div>

Hi,  
I want to build an MPO using a custom defined ITensor operator. The MPO can be constructed using  
MPO(::sites, ::list of op) (reference to definition : [https://itensor.github.io/ITensors.jl/dev/MPSandMPO.html](https://itensor.github.io/ITensors.jl/dev/MPSandMPO.html)).  
However, the caveat is that the list of operators are two-body. For example, for spin half systems, this two body operator could be Z \otimes Z or the one shown below.

```auto
using ITensors
using LinearAlgebra

sites = siteinds("Qubit", 6)

function ITensors.op(::OpName"TwoSiteOp", ::SiteType"Qubit"; θ1, θ2)
    
    M = cos(θ1)*Diagonal(ones(4)) - im*sin(θ2)*[1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1]

    return M
end

```

Now, the question is how can one construct an MPO for the defined ‘siteinds’, sites?  
The op list will be made of various operators with different \theta 1 and \theta 2.

Any suggestion for this will be very helpful.

---

<div class="post-metadata">

**Author:** ![ryanlevy](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ryanlevy/32/150_2.png) [@ryanlevy](https://itensor.discourse.group/u/ryanlevy)\
**Post date:** [August 13, 2024, 7:58pm UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/2 "2024-08-13T19:58:30Z")

</div>

See Miles’ answer from here:

> [@Defining two site local operators as part of a local Hamiltonian](https://itensor.discourse.group/t/defining-two-site-local-operators-as-part-of-a-local-hamiltonian/386/2):
>
> It’s an interesting thing to do, and would be a good feature for us to have, but unfortunately right now the AutoMPO system (since renamed the OpSum system) doesn’t support two-site operators like “CZ” currently. Internally what it would have to do is to compute the SVD of this operator to write it as a sum of products of single-site operators, and then put these into the sum of terms (like a rewriting rule). So it’s something we could add but would take a bit of development. As a workaround, y…

Essentially the OpSum system can’t do that at the moment, however if you can write it as a string (O\_i \otimes O\_j \dots ) or a sum of strings as in the post then that can be used. And for generic 2 site matrices you would have to decompose it by hand into that format

---

<div class="post-metadata">

**Author:** ![Bhargava\_Balaganchi](https://avatars.discourse-cdn.com/v4/letter/b/a5b964/32.png) [@Bhargava\_Balaganchi](https://itensor.discourse.group/u/Bhargava_Balaganchi)\
**Post date:** [August 14, 2024, 9:37am UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/3 "2024-08-14T09:37:06Z")

</div>

Hi Ryan,  
Thank you very much for the reply. To clarify, the two-site operator I mentioned above, e^{-i \theta\_1 \mathbb Z\_{1} \otimes \mathbb Z\_{2}}, can be decomposed as a product of single-site tensors as  
\begin{pmatrix} \cos(\theta\_1) \mathbb{I}\_{1} &-i\sin(\theta\_1) \mathbb Z\_{1} \end{pmatrix} \otimes \begin{pmatrix} \mathbb{I}\_{2} \\ \mathbb Z\_{2} \end{pmatrix}.  
Is there a way to construct an MPO by using such tensors as input?

---

<div class="post-metadata">

**Author:** ![ryanlevy](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ryanlevy/32/150_2.png) [@ryanlevy](https://itensor.discourse.group/u/ryanlevy)\
**Post date:** [August 14, 2024, 11:52am UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/4 "2024-08-14T11:52:52Z")

</div>

If you have the operator O = \cos(\theta\_1) I\_1 \otimes I\_2 -i\sin(\theta\_1) Z\_1 \otimes Z\_2 then the corresponding OpSum call would be

```julia
os = OpSum()
os += cos(θ1), "I",1,"I",2 # the second I is not needed really
os += -1im*sin(θ2),"Z",1,"Z",2
O = MPO(os,sites)

```

---

<div class="post-metadata">

**Author:** ![Bhargava\_Balaganchi](https://avatars.discourse-cdn.com/v4/letter/b/a5b964/32.png) [@Bhargava\_Balaganchi](https://itensor.discourse.group/u/Bhargava_Balaganchi)\
**Post date:** [August 14, 2024, 4:20pm UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/5 "2024-08-14T16:20:49Z")

</div>

Hi Ryan,  
thank you so much! There might still be a miscommunication though. Just to be sure: I need to build an MPO for the operator  
e^{-i\theta\mathbb Z\_{1} \otimes \mathbb Z\_{2}} \otimes e^{-i\theta\mathbb Z\_{3} \otimes \mathbb Z\_{4}} \otimes ....... \otimes e^{-i\theta\mathbb Z\_{N-1} \otimes \mathbb Z\_{N}} .  
This is a product rather than sum of two-site operators (unlike the link with Miles’ answer).  
Taking your suggestion to use OpSum, I tried the following code :

```auto
using ITensors

N = 8
sites = siteinds("Qubit", N)
U = MPO(sites)

θ = 0.5
os = OpSum()
os += cos(θ), "I", 1, "I", 2
os += -1im * sin(θ), "Z", 1, "Z", 2

for j in 1:2:N-1
    U[j:j+1] = MPO(os, sites[j:j+1])
end

```

However, this does not seem to produce the correct results. Could you help me correct it?

---

<div class="post-metadata">

**Author:** ![ryanlevy](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/ryanlevy/32/150_2.png) [@ryanlevy](https://itensor.discourse.group/u/ryanlevy)\
**Post date:** [August 14, 2024, 5:18pm UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/6 "2024-08-14T17:18:34Z")

</div>

You’ll need to decompose your operator into a list of single site operators  
 O = \sum\_n O\_1^n \otimes O\_2^n \otimes \dots, which for each n you would write the `os+=` line as above. You’d then call `MPO` just once to form the MPO.

Then one trick you can do, if you want say all "Z"s like A\*Z\_1 \otimes Z\_2 \otimes \dots \otimes Z\_N where A is a constant, then you can write

```julia
A = 0.1 # etc
term = [("Z",i) for i=1:N]
add!(os, A, (term...)...)

```

---

<div class="post-metadata">

**Author:** ![Bhargava\_Balaganchi](https://avatars.discourse-cdn.com/v4/letter/b/a5b964/32.png) [@Bhargava\_Balaganchi](https://itensor.discourse.group/u/Bhargava_Balaganchi)\
**Post date:** [August 15, 2024, 8:42am UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/7 "2024-08-15T08:42:09Z")

</div>

Thank you very much.

---

<div class="post-metadata">

**Author:** ![corbett5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/corbett5/32/886_2.png) [@corbett5](https://itensor.discourse.group/u/corbett5)\
**Post date:** [October 17, 2024, 9:12pm UTC](https://itensor.discourse.group/t/mpo-construction-from-two-site-custom-op/1924/8 "2024-10-17T21:12:27Z")

</div>

If you expand everything out in terms of single site operators as @ryanlevy suggested then you’ll have 2^{N / 2} different products which won’t scale well. I think you, @Bhargava_Balaganchi, had the right idea in your code example but didn’t set the elements of U correctly.

```julia
using ITensors

N = 8
sites = siteinds("Qubit", N)
U = MPO(sites, "I")

θ = 0.5
os = OpSum()
os += cos(θ), "I", 1, "I", 1
os += -1im * sin(θ), "Z", 1, "Z", 2

for j in 1:2:N-1
  Ui = MPO_new(os, sites[j:(j + 1)])
  U[j] = Ui[1]
  U[j + 1] = Ui[2]
end

# The purpose of orthogonalize! is to insert the missing link dims of size 1. 
orthogonalize!(U, 1)
@show linkdims(U)

```

```auto
linkdims(U) = [2, 1, 2, 1, 2, 1, 2]

```

Alternatively, a more general approach is to construct MPOs for the individual terms and then multiply them together. For this you’ll want to use the `ITensorMPOConstruction` library since `ITensors`’ implementation adds extra links for the starting and ending identity operator which can really add up when you multiply a bunch of MPOs together.

```julia
using ITensors
using ITensorMPOConstruction

N = 8
θ = 0.5

sites = siteinds("Qubit", N)
U = MPO(sites, "I")

for j in 1:2:N-1
  os = OpSum()
  os += cos(θ), "I", j, "I", j + 1
  os += -1im * sin(θ), "Z", j, "Z", j + 1
  Ui = MPO_new(os, sites)
  U = apply(Ui, U; alg="naive", truncate=false)
end

@show linkdims(U)

```

```auto
linkdims(U) = [2, 1, 2, 1, 2, 1, 2]

```

If you were to use `ITensors`’ MPO construction algorithm you would get

```auto
linkdims(U) = [32, 16, 32, 16, 32, 16, 32]

```

and in general find the the link dimension grows exponentially with N.