# Bug in multiplying MPOs and tracing MPOs

**URL:** <https://itensor.discourse.group/t/bug-in-multiplying-mpos-and-tracing-mpos/998>\
**Category:** ITensor Julia Questions\
**Tags:** julia, mpo\
**Created:** [July 3, 2023, 7:29pm UTC](https://itensor.discourse.group/t/bug-in-multiplying-mpos-and-tracing-mpos/998 "2023-07-03T19:29:41Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Brian-Khor](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/brian-khor/32/297_2.png) [@Brian-Khor](https://itensor.discourse.group/u/Brian-Khor)\
**Post date:** [July 3, 2023, 7:29pm UTC](https://itensor.discourse.group/t/bug-in-multiplying-mpos-and-tracing-mpos/998/1 "2023-07-03T19:29:41Z")

</div>

Hi admin,

I would like to implement a code that, given an MPS | \psi\_{MPS} \rangle, and assuming that I have correctly obtained density matrix \rho = | \psi \rangle \langle \psi | and the reduced density operator \rho\_A for some bipartition A/B (I use 60 sites below and use exactly half bipartition), calculate the quantity

Tr( \rho\_{A, Q}^2 ) = \int\_{- 2 \pi}^{2 \pi} d \lambda ( - \frac{| \lambda |}{(2 \pi)^2} + \frac{1}{2 \pi} ) Tr( \rho\_A e^{i \lambda Q} \rho\_A e^{- i \lambda Q} )

where Q = \sum\_{i=1}^{L\_A -1} \sigma\_i^z \sigma\_{i+1}^z should ideally be an MPO so that we can apply its exponential to the reduced density matrix \rho\_A effectively.

I have 2 questions:  
(1) I have error when I’m trying to apply MPOs e^{i \lambda \sigma\_i^z \sigma^z\_{i+1}} to \rho\_A. I think I cast both of them as MPOs but I suspect my code earlier might have error.

(2) For full trace, is there a built in function in iTensor Julia that takes in a final MPO and spits out the trace?

(3) Is my way, i.e., orthogonalize and then contract (for both reduced density operator and full trace) the most efficient way?

The code:

```julia
# (1) Construct Half Chain Reduced Density Matrix via contraction
    
# Orthogonalize so that we are in left-right canonical form
orthogonalize!(psi, 30)
    
# Density matrix as outer product of mixed canonical MPS
rho = outer(psi', psi)
sites = siteinds(rho)
    
L = ITensor(1.0)
for i=31:N
        # Need to check if delta is needed here
        L *= rho[i] * delta(sites[i][1], sites[i][2])
        
        # This printing check no of indices as sanity check
        # @show inds(L)
        # pause()
end
    
# Absorb the partial trace result to site 30
rho[30] *= L
            
# Cast the partial traced rho into MPO form
 rho_reduced = MPO(30)
 for i=1:30
     rho_reduced[i] = rho[i]
 end
    
# (2) Calculate Half Chain Symmetry Resolved Renyi Entropy
# via numerical integration
    
 Renyi_Sym = 0.0
 lambda = -2*pi
 prev_tr = 0.0
 curr_tr = 0.0
 new_rho_reduced = copy(rho_reduced)
    
# k = No of integration steps/splits
for k=1:101
    # Apply exponential of Q operator to rho_reduced for lambda
    # via rho_A e^(-i lambda Q) rho_A e^(i lambda Q)
            
    temp1_rho_reduced = deepcopy(rho_reduced)
    temp2_rho_reduced = deepcopy(rho_reduced)
        
    # Apply MPO of e^{i sigma^z_i sigma^z_{i+1}} for sites in A
    for alpha=1:29
        s1 = s[alpha]
        s2 = s[alpha+1]
        os = 4 * op("Sz", s1) * op("Sz", s2)
        szszp = exp(1.0im * lambda * os)
        szszm = exp(-1.0im * lambda * os)
        SzSzp = MPO([szszp])
        SzSzm = MPO([szszm])
            
        # Error here
        temp1_rho_reduced = apply(SzSzp, temp1_rho_reduced)
        temp2_rho_reduced = apply(SzSzm, temp2_rho_reduced)
    end
                    
    new_rho_reduced = apply(temp1_rho_reduced, temp2_rho_reduced)
    orthogonalize!(new_rho_reduced, 30)
    temp_sites = siteinds(new_rho_reduced)
            
    # Full trace for lambda
    for site=1:30
        new_rho_reduced[site] = new_rho_reduced[site] * delta(temp_sites[site][1], temp_sites[site][2])
    end
        
    curr_tr = (- (abs(lambda)/(2*pi)^2) + (1/(2*pi)))*(new_rho_reduced)

    # Perform trapezoid integration
    if k>1
        Renyi_Sym += 0.5*(prev_tr + curr_tr)*((2*pi)/50)
    end
            
    # Update
    lambda += (2*pi)/100
    prev_tr = curr_tr
end
    
Renyi_Sym = - log2(Renyi_Sym)

```

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [September 4, 2023, 2:46pm UTC](https://itensor.discourse.group/t/bug-in-multiplying-mpos-and-tracing-mpos/998/2 "2023-09-04T14:46:29Z")

</div>

Sorry for the very slow reply from us. (One tip also about asking questions on here is to just ask them one at a time.)

About question (2): yes, you can use the function `tr` to trace an MPO. It was not documented so I have just submitted a PR on Github to add documentation for it.

When you ask if (3) is the most efficient way, what objective are you referring to? (I.e. most efficient way to do which thing?)

---

<div class="post-metadata">

**Author:** ![Brian-Khor](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/brian-khor/32/297_2.png) [@Brian-Khor](https://itensor.discourse.group/u/Brian-Khor)\
**Post date:** [October 23, 2023, 7:30pm UTC](https://itensor.discourse.group/t/bug-in-multiplying-mpos-and-tracing-mpos/998/3 "2023-10-23T19:30:03Z")

</div>

Hi Miles, thanks for the reply. I managed to resolve the bug, after being careful with priming the outer products, contracts, etc. Now looking back I realised question (3) makes no sense.  
I ended up using tr function like you said, though I will also appreciate a documentation of the tr function.