Error in partial trace

[Rafael Barfknecht]

Hello all,

I have been trying to calculate the entropy for the reduced density matrix in a two mode system (basically a double well problem in a fixed number of excitations).

Since I did not find a way to directly calculate the entropy of reduced density matrices, I tried taking the partial trace of the full density matrix. When I have a single particle in the system, it seems to work fine since the system reproduces a two-qubit problem.

When I increase the number of particles (n=2 for instance), I get errors such as the following:

----------------------------------------

Traceback (most recent call last):

  File "/Users/rafaeleb/Desktop/question.py", line 64, in <module>

    rho1 = rhot.ptrace(0)

  File "/opt/anaconda3/envs/qutip-env/lib/python3.9/site-packages/qutip/qobj.py", line 1382, in ptrace

    return _ptrace_dense(self, sel)

  File "/opt/anaconda3/envs/qutip-env/lib/python3.9/site-packages/qutip/qobj.py", line 2234, in _ptrace_dense

    rhomat = np.trace(Q.full()

ValueError: cannot reshape array of size 64 into shape (3,3,3,3)

----------------------------------------

So my question is, can I somehow get the reduced density matrix in a system with a larger number of particles or is this function restricted to the 2-qubit system? See a sample of the code below.

Best regards,

Rafael

----------------------------------------

from qutip import *

import numpy as np

import matplotlib.pyplot as plt

import sys

nb = 2 #total number of bosons

nH = nb + 1 #dimension of each well

a = enr_destroy([nH,nH], excitations = nH) #list of annihilation operators acting on each site with fixed number of excitations

n1 = a[0].dag()*a[0]

n2 = a[1].dag()*a[1]

h1 = n1*(n1-1)

h2 = n2*(n2-1)

U0 = 0.1 #local interaction

t = 1. #tunneling rate

h0 = h1 + h2

h12 = a[0].dag()*a[1] + a[0]*a[1].dag()

H = - t*h12 + 0.5*U0*h0 #full Hamiltonian

e,v = H.eigenstates()

psi0 = enr_fock([nH,nH], excitations = nH, state = (nb,0))

tfinal = 100

steps = tfinal*10

times = np.linspace(0, tfinal, steps)

g1 = 0.

g2 = 0.1

#dissipation operators

lk1 = np.sqrt(g1)*a[0]

lk2 = np.sqrt(g1)*a[1]

#decoherence operators

dc1 = np.sqrt(g2)*n1

dc2 = np.sqrt(g2)*n2 

opts = Options()

opts.store_states = True

result = mesolve(H, psi0, times, [dc1,dc2], [n1,n2], options = opts, progress_bar = True)

fig, ax = plt.subplots()

n1exp = result.expect[0]

n2exp = result.expect[1]

dms = result.states

vnelistt = []

for j in range(len(dms)):

    print('Density matrix number '+str(j), end='\r')    

    rhot = dms[j]

    vnet = entropy_vn(rhot)

    vnelistt.append(vnet)

print(rhot)

rho1 = rhot.ptrace(0)

[Neill Lambert]

Hi Rafael,

The problem might be that the default ptrace (like many standard qutip functions)  doesn't work on ENR states.

There is a special ENR_ptrace function for enr states provided at the end of this example notebook

https://github.com/qutip/qutip-notebooks/blob/master/examples/excitation-number-restricted-states-jc-chain.ipynb

That should work, but it is not particularly optimized or tested, I just bodged it together after someone asked for it a few years back.    

Perhaps it might be another good thing to improve upon for the ENR stuff in qutip proper.

Thanks

Neill

--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/371e40e0-b7e6-45cf-bb9d-586e71cd69b9n%40googlegroups.com.

--

-----------------------------------
Senior Research Scientist
RIKEN, Japan
Research Homepage

[Rafael Barfknecht]

Hi Neil,

Thanks for the fast reply! I added your ENR_ptrace function to my code and it seems to work just fine! I will run more tests and write back if I find any issues. 

A very basic question: I would like to add this function to the library to use like the other enr functions. How do I go about doing that? I installed QuTip via anaconda, but I cannot find the files where the functions are defined.

Best,

Rafael

[Boxi Li]

Hi Rafael,

The other enf functions are in operators.py and states.py. 

If you want to add this to the qutip library and use it as part of qutip, you need to hack the installed package: locate your conda path and edit the files there. However, if you want to do things properly, you will need to install qutip in the developer mode. This means that you will need to download the source code on our GitHub add install it from the source (see http://qutip.org/docs/latest/installation.html#direct-setuptools-source-builds and use the command `python setup.py develop`). There might be a few detours and you will probably need C++ compiler to make everything work.

We are also happy if you are willing to make this part of qutip and everyone else can use it. We have a written guide on contributing: http://qutip.org/docs/latest/development/contributing.html, if you are interested.

Best

Boxi

[Rafael Barfknecht]

Hi Boxi,

Thank you for the reply. I did find the files where I can add the functions. For now I will use this solution, but I would be happy to contribute where I can in the future. What I can say at this point is that Neil's function seems to work fine for my purposes, so I believe it would make sense for it to be included along with the other enr functions in a future QuTip update.

Best wishes,

Rafael